Graphical user interfaces for the R language are easy to use and getting more powerful all the time. Here is my updated comparison of jamovi, JASP, BlueSky Statistics (free & Pro), Rattle, RKWard, R-Instat, R AnalyticFlow, and R Commander.

With so many detailed reviews of Graphical User Interfaces (GUIs) for R available, which should you choose? It’s not too difficult to rate them based on the number of features they offer, so I’ll start there. Then, I’ll follow with a brief overview of each.

I’m basing the counts on the number of dialog boxes in each category of the following categories:

• Ease of Use
• General Usability
• Graphics
• Analytics
• Reproducibility

This data is trickier to collect than you might think. Some software has fewer menu choices, depending instead on more detailed dialog boxes. Studying every menu and dialog box is very time-consuming, but that is what I’ve tried to do to keep this comparison trustworthy.

Perhaps the biggest flaw in this methodology is that every feature adds only one point to each GUI’s total score. I encourage you to download the full dataset and consider which features are most important to you. If you decide to make your own graphs with a different weighting system, I’d love to hear from you in the comments below. Thomas Langkamp did an excellent comparison focused on teaching statistics. You can read it here.

Ease of Use

For ease of use, I’ve defined it primarily by how well each GUI meets its primary goal: avoiding code. They get one point for each of the following abilities, which include being able to install, start, and use the GUI to its maximum effect, including publication-quality output, without knowing anything about the R language itself. Figure 1 shows the result. R Commander is abbreviated Rcmdr, and R AnalyticFlow is abbreviated RAF. The commercial BlueSky Pro comes out on top by a slim margin, followed closely by JASP and RKWard. None of the GUIs achieved the highest possible score of 15, so there is room for improvement.

The ease of use features are:

• Installs without the use of R
• Starts without the use of R
• Remembers recent files
• Hides R code by default
• Uses its full capability without needing R code
• Data editor included
• Pub-quality tables w/out R code steps
• Simple menus that grow as needed
• Table of Contents to ease navigation
• Variable labels ease identification in the output
• Easy to move blocks of output
• Ease reading columns by freezing headers of long tables
• Accepts data pasted from the clipboard
• Easy to move header row of pasted data into the variable name field
• Includes AI interface to ease interpretation of output

General Usability

This category is dominated by data-wrangling capabilities, which data scientists and statisticians spend most of their time on. It also includes various types of data input and output. We see in Figure 2 that both BlueSky versions and R-Instat come out on top, not just due to their excellent selection of data-wrangling features but also for their use of the rio package for importing and exporting files. The rio package combines the import/export capabilities of many other packages, and it is easy to use. I expect the other GUIs will eventually adopt it, raising their scores by around 20 points.

The general usability features are:

• Operating systems (how many)
• Import data file types (how many)
• Import from databases (how many)
• Export data file types (how many)
• Languages displayable in UI (how many, besides English)
• Easy to repeat any step by groups (split-file)
• Multiple data files open at once
• Multiple output windows
• Multiple code windows
• Variable metadata view
• Variable types (how many)
• Variable search/filter in dialogs
• Variable sort by name
• Variable sort by type
• Variable move manually
• Model Builder (how many effect types)
• Magnify GUI for teaching
• R code editor
• Comment/uncomment blocks of code
• Package management (comes with R and all packages)
• Output: word processing features
• Output: R Markdown
• Output: LaTeX
• Data wrangling (how many)
• Transform across many variables at once (e.g., row mean)
• Transform down many variables at once (e.g., log, sqrt)
• Assign factor labels across many variables at once
• Project saves/loads data, dialogs, and notes in one file

Graphics

This category consists mainly of the number of plot types each software offers. However, the other items can be very important to completing your work. They should add more than one point to the graphics score, but I scored them one point since some will view them as very important while others might not need them at all. Be sure to see the full reviews or download the Excel file if those features are important to you. Figure 3 shows the total graphics score for each GUI. jamovi has a solid lead in this category. R-Instat’s ability to layer any geometric object, or “geom,” on top of another could substantially increase its count, but I don’t include that since it requires knowledge of R code.

When studying all the graphs in this article, it’s important to consider the difference between the relative and absolute performance. For example, relatively speaking, R Commander is not doing well here, but it does offer over 25 types of plots! That might be fine for your needs.

The graphics features are:

• Types of graphs (how many)
• Small multiples (i.e., facets)
• Large multiples (i.e., by group / split-file)
• Export graphics formats (how many)
• Interactive Plot Builder
• Interactive Plot Editor (e.g., add annotations, arrows, etc.)
• Able to combine any graphs into a single image

Continued…

The PROJ format for defining Coordinate Reference Systems (CRS) is becoming deprecated, as PROJ strings cannot accurately store modern datum transformations, coordinate epoch information, or complex axis order rules. However, that format is much easier and more practical for humans to type and read, and it is still often used as a legacy convenience. I, for one, use it often for centering a map projection on the mapped region, thus avoiding the distortion that occurs away from the projection center. For example, imagine I want to map New Zealand in a Lambert Azimuthal Equal-Area projection in R:

library(terra)library(geodata)nz <- gadm(country = "New Zealand", level = 0, path = tempdir())nz_prj <- project(nz, "+proj=laea")plot(nz_prj)

The map is heavily distorted, because it is far away from the center of the projection, which is 0,0 by default. With the PROJ format, we can easily set the center wherever we want, e.g.:

nz_centr <- crds(centroids(aggregate(nz)))nz_centr  # 172.6104 -41.82592nz_prj <- project(nz, "+proj=laea +lat_0=-41.8 +lon_0=172.6")# (this can be further automated to avoid the hard-coding,# but I keep it like this here for simplicity) plot(nz_prj)

Much better! We now have an equal-area map of New Zealand without visibly distorting the shape, and without having to look for specific equal-area CRSs normally used over there (and again every time we need to map another region of the world). This works as long as we use the WGS84 ellipsoid rather than a local datum, as those are currently not supported when using PROJ. So, I’m still using this regularly, and apprehensive about PROJ being deprecated… If anyone knows of a non-PROJ way of easily defining user-centered projections like this, please let me know!

After an extended beta period, JAGS 5.0.0 is now available. The release is more of a process than a single event, but we have started by uploading the source tarball and the Windows installer to SourceForge. The macOS installer will be along hopefully next week.

Thanks to Jeroen Ooms, Windows binary installers are available from the r-windows project on GitHub. This includes binaries for aarch64/Arm64.

Changes since the beta release

• Better multi-core support on macOS (More on this in a future post).
• Fixed a delayed rejection bug in the MNormal sampler.
• Added a new elliptical slice sampler contributed by Marcel Jonker.
• The dic.samples function could throw a logic error (See ticket #66). I am unable to reproduce this but JAGS 5.0.0 includes a workaround that should avoid triggering this bug.
• Fixed a reproducibility bug in BaseRNGFactory.

I would not normally add a new sampler at this late stage, but the elliptical slice sampler was ready to go and performs surprisingly well for a non-adaptive sampler. If you experience any trouble with this new sampler then you can turn it of with

set.factory("bugs::EllipticalSlice", type="sampler", state=FALSE)

JAGS will then fall back on the adaptive random walk Metropolis-Hastings sampler MNormal.

The rjags package and CRAN migration

The current rjags package on CRAN (version 4-17) will not work with JAGS 5.0.0. You need version 5-x but for the time being this will only be available from SourceForge.

• Source tarball
• Windows binary (x86_64)
• macOS binary

The reason for not uploading the rjags package to CRAN is that it currently breaks a lot of reverse dependencies. I have been quite proactive in contacting the maintainers of these package and providing patches, but we still need to give them more time.

We are monitoring the reverse dependency situation and will be liaising with the CRAN maintainers.

Preview of upcoming posts

In order to help users and developers with the transition to JAGS version 5, a series of “What’s new in JAGS 5.0.0” blog posts are planned covering the following topics.

• Initial values
• Reproducibility
• Multicore support
• ABI compatibility
• New diagnostics

These will include guidance for package maintainers to make the transition to JAGS 5.0.0.

I’ve put together a short guide: How Do I Build a Model? A starting guide for new modellers. It’s free, online, and aimed squarely at the graduate student who has a research question and is staring at a blank script wondering where to even begin.

You can read it online at How Do I Build a Model?, or download a pdf.

I’ve never settled into one corner of ecological modelling. I’m not the Bayesian guy, or the foodweb model guy, or the machine learning guy, though I’ve dabbled in all three. That haphazard career has given me a broad view across a lot of different modelling traditions, and I keep noticing the same meta-principles showing up regardless of the discipline. This guide is my attempt to write those down.

It’s not a textbook on any particular method. Once you know what type of model you’re building, go read a book by a discipline expert (I’ve listed some favourites in the final chapter). This guide is for the step before that: turning a question into a model in the first place.

In their paper “A Synchronized Multi‑IMU Wearable System for Tracking of Joint‑Angles in Sports Motion Analysis” (arXiv:2607.26027v1), Samarasekera and colleagues set out to solve a very practical problem: how to reliably measure joint angles in dynamic sports movements using wearable IMUs. Their goal was to design a synchronized pipeline that could filter out noise, correct drift, and normalize ranges so that biomechanical signals could be trusted in real‑world conditions.

Our goal is different but conceptually parallel: to apply the same pipeline logic to financial time‑series data. Instead of knee flexion angles, we analyze the comparative dynamics of Gold and Bitcoin. The challenge is similar — financial returns are noisy, prone to drift, and difficult to compare directly. By adapting their pipeline, we aim to extract latent trends, mitigate drift, normalize ranges, and highlight relative divergence between assets.

The pipeline consists of four sequential stages:

1. Indirect Kalman Filtering (IKF) – extracting latent trends from noisy log returns.
   • State equation:

$$x_t=F\cdot x_{t-1}+w_t,w_t\sim N(0,Q)$$

• Observation equation:

$$y_t=H\cdot x_t+v_t,v_t\sim N(0,R)$$

• IKF separates signal from noise, with $Q$ and $R$ learned via maximum likelihood.

2. Drift Mitigation via High‑Pass Filtering – removing low‑frequency bias.
   • High‑pass operator:

$$x_t^{HP}=x_t-\frac{1}{N}\sum _{i=1}^N{x}_{t-i}$$

• This ensures long‑term drift does not dominate short‑term dynamics.

3. Range Normalization – mapping filtered values into a comparable scale.
   • Normalization function:

$$z_t=\frac{x_t^{HP}-\min (x^{HP})}{\max (x^{HP})-\min (x^{HP})}\cdot 100$$

• This allows Gold and Bitcoin to be compared on the same 0–100 scale.

4. Relative Notation (Directional Divergence) – highlighting when assets move in opposite directions.
   • Divergence indicator:

$$D_t=\text{sign}(r_t^{Gold})\ne \text{sign}(r_t^{BTC})$$

• If true, shaded regions show divergence; if false, alignment.

Why This Matters

• For the biomechanics researchers, the pipeline meant trustworthy joint‑angle tracking.
• For us, the same architecture means trustworthy comparative trend analysis between financial assets.
• In both cases, the pipeline’s strength lies in its modularity: IKF for latent signal extraction, high‑pass filtering for drift control, normalization for comparability, and relative notation for divergence detection.

✦ R Packages in Our Pipeline

To make this pipeline reproducible and transparent, we rely on a set of R packages. Each package plays a specific role in the architecture:

1. tidyverse
   • Purpose: Data manipulation and wrangling.
   • Functions: mutate, rename, inner_join, drop_na.
   • Why: Provides a clean grammar for transforming raw financial data into structured time‑series.
2. tidyquant
   • Purpose: Financial data acquisition and transformation.
   • Functions: tq_get (download asset prices), tq_transmute (compute log returns).
   • Why: Bridges tidyverse with financial APIs, enabling reproducible market data pipelines.
3. KFAS
   • Purpose: State‑space modeling and Kalman filtering.
   • Functions: SSModel (define state‑space model), fitSSM (MLE parameter estimation), KFS (Kalman smoothing).
   • Why: Implements the mathematical backbone of IKF, separating latent trends from noise.
4. zoo
   • Purpose: Rolling window operations.
   • Functions: rollmean.
   • Why: Used for high‑pass filtering by subtracting rolling averages, mitigating drift.
5. scales
   • Purpose: Normalization and rescaling.
   • Functions: rescale.
   • Why: Maps filtered trends into a 0–100 range, enabling comparability across assets.
6. ggbraid
   • Purpose: Divergence visualization.
   • Functions: geom_braid.
   • Why: Highlights periods of divergence between Gold and Bitcoin with shaded regions.
7. ggplot2 (via tidyverse)
   • Purpose: Visualization.
   • Functions: geom_line, labs, theme_classic.
   • Why: Provides the plotting framework for presenting normalized high‑pass Kalman trends.

library(tidyverse)
library(tidyquant)
library(KFAS)
library(ggbraid)

# --- Step 1: Log returns ---
gold <- tq_get("GC=F", from = "2024-01-01") %>%
  tq_transmute(select = close,
               mutate_fun = periodReturn,
               period = "daily",
               type = "log")

btc <- tq_get("BTC-USD", from = "2024-01-01") %>%
  tq_transmute(select = close,
               mutate_fun = periodReturn,
               period = "daily",
               type = "log")

df <- gold %>%
  rename(ret_gold = daily.returns) %>%
  inner_join(btc %>% rename(ret_btc = daily.returns), by="date") %>%
  drop_na()

# --- Step 2: IKF (Kalman filter) ---
model_gold <- SSModel(df$ret_gold ~ SSMtrend(1, Q = list(NA)), H = NA)
model_btc  <- SSModel(df$ret_btc  ~ SSMtrend(1, Q = list(NA)), H = NA)

fit_gold <- fitSSM(model_gold, inits = c(log(var(df$ret_gold)), log(var(df$ret_gold))))
fit_btc  <- fitSSM(model_btc,  inits = c(log(var(df$ret_btc)),  log(var(df$ret_btc))))

kf_gold <- KFS(fit_gold$model, smoothing = c("state"))
kf_btc  <- KFS(fit_btc$model,  smoothing = c("state"))

trend_gold <- as.numeric(kf_gold$alphahat)
trend_btc  <- as.numeric(kf_btc$alphahat)

# --- Step 3: High-pass filter (remove low-frequency drift) ---
hp_gold <- trend_gold - zoo::rollmean(trend_gold, k=30, fill=NA, align="right")
hp_btc  <- trend_btc  - zoo::rollmean(trend_btc,  k=30, fill=NA, align="right")

# --- Step 4: Range normalization ---
gold_norm <- scales::rescale(hp_gold, to=c(0,100))
btc_norm  <- scales::rescale(hp_btc,  to=c(0,100))

# --- Step 5: Relative notation (directional divergence) ---
df_norm <- df %>%
  mutate(
    gold_norm = gold_norm,
    btc_norm  = btc_norm,
    sign_gold = sign(ret_gold),
    sign_btc  = sign(ret_btc),
    divergence = sign_gold != sign_btc
  )

# --- Plot ---
ggplot(df_norm, aes(x=date)) +
  geom_line(aes(y=gold_norm, color="Gold"), linewidth=1.2) +
  geom_line(aes(y=btc_norm, color="Bitcoin"), linewidth=1.2) +
  geom_braid(aes(ymin=gold_norm, ymax=btc_norm, fill=divergence), alpha=0.3) +
  scale_color_manual(values=c("Gold"="#FFD700", "Bitcoin"="#1E3A8A"),
                     name = "") +
  scale_fill_manual(values=c("TRUE"="#FF4C4C", "FALSE"="#87CEFA"),
                    name = "Divergence") +
  scale_y_continuous(limits=c(0,100)) +
  labs(title="Gold vs Bitcoin: Kalman-Filtered High-Pass Trends",
       subtitle="Directional divergence highlighted with shaded regions",
       x="", y="Normalized High-Pass Kalman Trend") +
  theme_classic(base_family = "Roboto Slab") +
  theme(
    legend.position="bottom",
    legend.title=element_text(size=12, face="bold"),
    legend.text=element_text(size=10),
    plot.title=element_text(size=16, face="bold"),
    plot.subtitle=element_text(size=12),
    axis.title=element_text(size=12, face="bold"),
    axis.text=element_text(size=10),
    panel.grid.minor=element_blank()
  )

Conclusion

By adapting the IKF pipeline from the sports biomechanics paper, we demonstrate how state‑space filtering, drift mitigation, and normalization can reveal comparative dynamics among financial assets. This modular approach is reproducible and extensible to other domains where latent trends and divergence matter.

This is a post based on two other posts I wrote for LinkedIn. It has become quite long.

Run cell.
Inspect output.
Modify.
Run again.
Oops.
Restart kernel.
Run all.

for ...
    ...

We have all seen colourful signposts with great cities and their distances – how far each is from this current standpoint.

And fundamental question is, how correct these distances are? Or even, one should ask, how honest a particular signpost is, regarding the current standpoint.

So let’s put this utterly most fundamental question to the test And since it is vacation time, you can grab the R code and check it yourself

Image stumbling upon the post and wanted to check if this post is really relevant for this location. It can be a photo from other side of the world, hanging in the lobby of your holiday hotel And this sure ain’t honest signpost

The R code iterates through each city, finds matching points (where all circles meet) and searches for the nearest location. It uses OpenStreetMap (short: OSM) to extract all the geolocations and looks for possible intersections of all points

geocode_city <- function(city_name) {
 
  url <- modify_url(
    "https://nominatim.openstreetmap.org/search",
    query = list(
      q       = city_name,
      format  = "json",
      limit   = 1
    )
  )
  
  resp <- tryCatch(
    GET(url, user_agent("PointingSignFinder/1.0 (R script)")),
     error = function(e) {
       cat(" Error on endpoint")
       return(NULL)
     }
  )
  
  if (is.null(resp) || http_error(resp)) {
    cat(" Error on http\n")
    return(NULL)
  }
  
  result <- fromJSON(content(resp, as = "text", encoding = "UTF-8"))
  
  if (length(result) == 0) {
    cat("Error on result")
    return(NULL)
  }
  
  lat <- as.numeric(result$lat[1])
  lon <- as.numeric(result$lon[1])
  cat(sprintf(" found: %.4f°, %.4f°\n", lat, lon))
  
  Sys.sleep(1.1)
  
  list(lat = lat, lon = lon, display_name = result$display_name[1])
}


sign_location_finder <- function(cities,
                                 distances,
                                 tolerance = 50,
                                 coarse_res = 0.2,
                                 fine_res   = 0.02,
                                 nearby_radius = 200,
                                 nearby_min_pop = 100000) {
  
  stopifnot(length(cities) == length(distances))
  stopifnot(length(cities) >= 2)
  
  cat("Step 1: Geocoding cities\n")
  coords <- lapply(cities, geocode_city)
  
  failed <- which(sapply(coords, is.null))
  if (length(failed) > 0) {
    stop(sprintf("Location not founc: %s", paste(cities[failed], collapse = ", ")))
  }
  
  sign_data <- data.frame(
    city    = cities,
    lat     = sapply(coords, `[[`, "lat"),
    lon     = sapply(coords, `[[`, "lon"),
    dist_km = distances
  )


print(sign_data[, c("city", "lat", "lon", "dist_km")])
  
max_dist_deg <- max(sign_data$dist_km) / 111
lat_min <- min(sign_data$lat) - max_dist_deg - 5
lat_max <- max(sign_data$lat) + max_dist_deg + 5
lon_min <- min(sign_data$lon) - max_dist_deg * 2 - 5
lon_max <- max(sign_data$lon) + max_dist_deg * 2 + 5
  
 
lat_min <- max(lat_min, -85)
lat_max <- min(lat_max,  85)
lon_min <- max(lon_min, -180)
lon_max <- min(lon_max,  180)
  
cat(sprintf("\nStep 2: Search bounding box: lat [%.1f, %.1f], lon [%.1f, %.1f]\n",
              lat_min, lat_max, lon_min, lon_max))
  
  
  score_point <- function(plat, plon) {
    diffs <- sapply(seq_len(nrow(sign_data)), function(i) {
      d <- distHaversine(
        c(plon, plat),
        c(sign_data$lon[i], sign_data$lat[i])
      ) / 1000
      abs(d - sign_data$dist_km[i])
    })
    max(diffs)
  }
  
  
  cat(sprintf("\nStep 3: Search the grid (%.2f° resolution)...\n", coarse_res))
  grid <- expand.grid(
    lat = seq(lat_min, lat_max, by = coarse_res),
    lon = seq(lon_min, lon_max, by = coarse_res)
  )
  cat(sprintf("Checking %d grid points...\n", nrow(grid)))
  grid$score <- mapply(score_point, grid$lat, grid$lon)
  
  candidates <- grid[grid$score <= tolerance, ]
  cat(sprintf(" Found %d candidate cells within ±%d km tolerance.\n",
              nrow(candidates), tolerance))
  
  coarse_best <- grid[which.min(grid$score), ]
  
  fine_grid <- expand.grid(
    lat = seq(coarse_best$lat - 1, coarse_best$lat + 1, by = fine_res),
    lon = seq(coarse_best$lon - 1, coarse_best$lon + 1, by = fine_res)
  )
  fine_grid$score <- mapply(score_point, fine_grid$lat, fine_grid$lon)
  best <- fine_grid[which.min(fine_grid$score), ]
  
  cat(sprintf("\n>>> Estimated  locations:\n"))
  cat(sprintf("    Latitude  : %.4f°\n", best$lat))
  cat(sprintf("    Longitude : %.4f°\n", best$lon))
  cat(sprintf("    Max error : ±%.1f km\n", best$score))
  cat(sprintf("    Google Maps: https://www.google.com/maps?q=%.4f,%.4f\n", best$lat, best$lon))
  
  cat(sprintf("\nStep 5: Nearby cities (within %d km, pop > %s)...\n",
              nearby_radius, format(nearby_min_pop, big.mark = ",")))
  
  world_cities <- world.cities
  nearby <- world_cities |>
    filter(pop > nearby_min_pop) |>
    mutate(
      dist_to_sign = distHaversine(
        cbind(long, lat),
        c(best$lon, best$lat)
      ) / 1000
    ) |>
    filter(dist_to_sign <= nearby_radius) |>
    arrange(dist_to_sign) |>
    select(name, country.etc, lat, long, pop, dist_to_sign) |>
    head(10)
  
  if (nrow(nearby) > 0) {
    cat("\n  Cities found:\n")
    print(nearby, digits = 4)
    cat(sprintf("\n  Nearest: %s, %s (%.1f km away)\n",
                nearby$name[1], nearby$country.etc[1], nearby$dist_to_sign[1]))
  } else {
    cat("No major locations found in vicinity!")
  }
  

}

So now, that we understand the useless problem, we can put this to the test

# sample 2 - real with correct distances calculater from#  https://www.distancefromto.net/  "air distance"# result must be Ljubljana, Sloveniaresult <- sign_location_finder(  cities    = c("Koper", "Celje", "Maribor", "Kranj"),  distances = c(83, 61, 104, 24),  tolerance = 20) 

In this case I have entered the correct air distances and there should be a location present – which is capital of Slovenia – Ljubljana.

I have inserted the real values of air distances between Ljubljana all all four cities: Koper, Celje, Maribor, Kranj using air distance calculator.

And this where the circles should meet

So if you are enjoying your vacation and you stumble upon a sign, check it and let me know, if it corresponds to your location of not

As always, the complete code is available on GitHub in  Useless_R_function repository and the file Honest_geographical_location_signposts.R is here.

Check the repository for future updates!

Stay healthy, hydrated and happy R-coding!

The package garchf wraps the extensive GARCH modeling machinery in rugarch behind a forecast-style interface, so that GARCH-based forecasts can be produced, plotted, and evaluated with the same conventions used throughout the forecast ecosystem (in the spirit of functions like forecast::thetaf()). The workhorse function, xgarchf(), fits a GARCH model to a univariate time series – with a choice of variance model (sGARCH, eGARCH, gjrGARCH, apARCH, iGARCH, …), ARMA order for the conditional mean, GARCH order, and conditional distribution – and returns a standard forecast object containing point forecasts, prediction intervals, and simulated paths. Because the output integrates with forecast’s S3 methods, plot() and other generics work out of the box, and the simulated paths (sims) make it straightforward to compute richer probabilistic forecasting metrics (CRPS, pinball loss, Winkler scores, etc.) via cross-validation. Below, we illustrate xgarchf() on Google’s daily closing stock price returns, and then run a small cross-validation exercise comparing a few GARCH variants using the crossvalidation package.

library(garchf)

Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo 

1 – Graphs

y <- diff(log(fpp2::goog200))

fit1 <- xgarchf(y, h = 20, model = "eGARCH")
fit2 <- xgarchf(y, h = 20, model = "sGARCH")
fit3 <- xgarchf(y, h = 20, model = "gjrGARCH")
fit4 <- xgarchf(y, h = 20, model = "iGARCH")

print(summary(fitted(fit1)))

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
0.001504 0.001504 0.001504 0.001504 0.001504 0.001504     
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
0.002783 0.002783 0.002783 0.002783 0.002783 0.002783     
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
0.001233 0.001233 0.001233 0.001233 0.001233 0.001233     
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
0.001485 0.001485 0.001485 0.001485 0.001485 0.001485 

    Ljung-Box test

data:  Residuals
Q* = 10.341, df = 10, p-value = 0.4111

Model df: 0.   Total lags used: 10

    Ljung-Box test

data:  Residuals
Q* = 10.341, df = 10, p-value = 0.4111

Model df: 0.   Total lags used: 10

    Ljung-Box test

data:  Residuals
Q* = 10.341, df = 10, p-value = 0.4111

Model df: 0.   Total lags used: 10

    Ljung-Box test

data:  Residuals
Q* = 10.341, df = 10, p-value = 0.4111

Model df: 0.   Total lags used: 10

2 - Cross-validation with the package crossvalidation

For crossvalidation, see https://github.com/Techtonique/crossvalidation.

spl_m5 <- function(predicted,
                   observed,
                   probs = c(0.005, 0.025, 0.165, 0.25,
                             0.5,
                             0.75, 0.835, 0.975, 0.995))
{
  if (is.null(predicted$sims)) {
    stop("predicted$sims is required")
  }

  # training series from forecast object
  train <- as.numeric(predicted$x)

  # M5 scaling denominator
  scale_denom <- mean(abs(diff(train)), na.rm = TRUE)

  if (scale_denom <= 0 || !is.finite(scale_denom)) {
    scale_denom <- 1
  }

  sims <- predicted$sims

  # quantile forecasts: rows = horizons, cols = probs
  qhat <- sapply(
    probs,
    function(p) apply(sims, 1, stats::quantile, probs = p, na.rm = TRUE)
  )

  observed <- as.numeric(observed)

  # pinball loss
  pinball_vec <- numeric(length(probs))

  for (j in seq_along(probs)) {

    u <- probs[j]
    q <- qhat[, j]

    pinball_vec[j] <- mean(
      ifelse(
        observed >= q,
        u * (observed - q),
        (1 - u) * (q - observed)
      ),
      na.rm = TRUE
    )
  }

  # scaled pinball loss for each quantile
  spl_vec <- pinball_vec / scale_denom

  # average scaled pinball loss
  mean_spl <- mean(spl_vec)

  out <- c(
    mean_SPL = mean_spl
  )

  names(spl_vec) <- paste0("SPL_", probs)

  c(out, spl_vec)
}

eval_metric95 <- function(predicted, observed)
{
  error <- observed - predicted$mean

  rmse <- sqrt(mean(error^2))
  mae <- mean(abs(error))

  # Only one interval returned
  lower <- predicted$lower
  upper <- predicted$upper

  coverage <- mean(
    observed >= lower & observed <= upper
  )

  alpha <- 0.05

  winkler <- ifelse(
    observed < lower,
    (upper - lower) + (2 / alpha) * (lower - observed),
    ifelse(
      observed > upper,
      (upper - lower) + (2 / alpha) * (observed - upper),
      (upper - lower)
    )
  )

  spl <- spl_m5(predicted, observed)

  c(
    RMSE = rmse,
    MAE = mae,
    Coverage95 = coverage,
    Winkler95 = mean(winkler),
    CRPS = mean(scoringRules::crps_sample(observed, predicted$sims)),
    mean_SPL = as.numeric(spl["mean_SPL"])
  )
}

eval_metric80 <- function(predicted, observed)
{
  error <- observed - predicted$mean

  me <- mean(error)
  rmse <- sqrt(mean(error^2))
  mae <- mean(abs(error))

  # Only one interval returned
  lower <- predicted$lower
  upper <- predicted$upper

  coverage <- mean(
    observed >= lower & observed <= upper
  )

  alpha <- 0.2

  winkler <- ifelse(
    observed < lower,
    (upper - lower) + (2 / alpha) * (lower - observed),
    ifelse(
      observed > upper,
      (upper - lower) + (2 / alpha) * (observed - upper),
      (upper - lower)
    )
  )

  spl <- spl_m5(predicted, observed)

  c(
    ME = me,
    RMSE = rmse,
    MAE = mae,
    Coverage80 = coverage,
    Winkler80 = mean(winkler),
    CRPS = mean(scoringRules::crps_sample(observed, predicted$sims)),
    mean_SPL = as.numeric(spl["mean_SPL"])
  )
}

res <- crossvalidation::crossval_ts(
  y = y,
  initial_window = 150,
  horizon = 10,
  fixed_window = FALSE,
  fcast_func = garchf::xgarchf,
  eval_metric = eval_metric95,
  fit_params = list(arma_order=c(0, 0), 
                    garch_order = c(1, 1),
                    model = "eGARCH",
                    level=95),
  show_progress = FALSE
)
print(summary(res))

res <- crossvalidation::crossval_ts(
  y = y,
  initial_window = 150,
  horizon = 10,
  fixed_window = FALSE,
  fcast_func = garchf::xgarchf,
  eval_metric = eval_metric95,
  fit_params = list(arma_order=c(0, 0), 
                    garch_order = c(1, 1),
                    model = "sGARCH",
                    level=95),
  show_progress = FALSE
)
print(summary(res))

res <- crossvalidation::crossval_ts(
  y = y,
  initial_window = 150,
  horizon = 10,
  fixed_window = FALSE,
  fcast_func = garchf::xgarchf,
  eval_metric = eval_metric95,
  fit_params = list(arma_order=c(0, 0), 
                    garch_order = c(1, 1),
                    model = "iGARCH",
                    level=95),
  show_progress = FALSE
)
print(summary(res))

      RMSE               MAE             Coverage95       Winkler95      
 Min.   :0.006459   Min.   :0.004896   Min.   :0.7000   Min.   :0.03441  
 1st Qu.:0.008001   1st Qu.:0.006266   1st Qu.:0.9000   1st Qu.:0.06153  
 Median :0.009135   Median :0.007061   Median :1.0000   Median :0.06466  
 Mean   :0.017189   Mean   :0.010613   Mean   :0.9475   Mean   :0.20090  
 3rd Qu.:0.022415   3rd Qu.:0.015066   3rd Qu.:1.0000   3rd Qu.:0.47706  
 Max.   :0.042741   Max.   :0.023054   Max.   :1.0000   Max.   :0.96625  
      CRPS             mean_SPL     
 Min.   :0.004154   Min.   :0.1170  
 1st Qu.:0.004922   1st Qu.:0.1362  
 Median :0.005311   Median :0.1481  
 Mean   :0.010039   Mean   :0.3380  
 3rd Qu.:0.017800   3rd Qu.:0.6763  
 Max.   :0.051416   Max.   :1.6607        

      RMSE               MAE             Coverage95       Winkler95      
 Min.   :0.006554   Min.   :0.005380   Min.   :0.7000   Min.   :0.03856  
 1st Qu.:0.008424   1st Qu.:0.006561   1st Qu.:0.9750   1st Qu.:0.06261  
 Median :0.009314   Median :0.007401   Median :1.0000   Median :0.07277  
 Mean   :0.017155   Mean   :0.010741   Mean   :0.9525   Mean   :0.17716  
 3rd Qu.:0.020081   3rd Qu.:0.013012   3rd Qu.:1.0000   3rd Qu.:0.34988  
 Max.   :0.042568   Max.   :0.022999   Max.   :1.0000   Max.   :0.49756  
      CRPS             mean_SPL     
 Min.   :0.004634   Min.   :0.1349  
 1st Qu.:0.005383   1st Qu.:0.1516  
 Median :0.005695   Median :0.1606  
 Mean   :0.008842   Mean   :0.2986  
 3rd Qu.:0.012903   3rd Qu.:0.4725  
 Max.   :0.019442   Max.   :0.7436        

    RMSE               MAE             Coverage95       Winkler95      
 Min.   :0.006556   Min.   :0.005514   Min.   :0.7000   Min.   :0.03580  
 1st Qu.:0.008306   1st Qu.:0.006239   1st Qu.:0.9750   1st Qu.:0.06307  
 Median :0.009119   Median :0.007007   Median :1.0000   Median :0.06478  
 Mean   :0.017029   Mean   :0.010484   Mean   :0.9525   Mean   :0.16725  
 3rd Qu.:0.018881   3rd Qu.:0.011711   3rd Qu.:1.0000   3rd Qu.:0.16760  
 Max.   :0.042713   Max.   :0.023064   Max.   :1.0000   Max.   :0.51546  
      CRPS             mean_SPL     
 Min.   :0.004851   Min.   :0.1363  
 1st Qu.:0.005317   1st Qu.:0.1510  
 Median :0.005700   Median :0.1591  
 Mean   :0.008720   Mean   :0.2928  
 3rd Qu.:0.009316   3rd Qu.:0.3126  
 Max.   :0.019649   Max.   :0.7586  

LLMs can produce a convincing media recommendation in a few seconds. The more useful question is whether the recommendation is correct: are the reach calculations right, are the assumptions visible, and does the plan fit the brief?

I work at Havas Media and I am the author of the Havas AI Media Quality Index (HAI-Q), a benchmark for media planning in the German market. I built it because generic knowledge tests are not enough for professional AI applications. A model can write a plausible paragraph about media strategy and still get a budget allocation, a target group definition, or a GRP calculation wrong.

This post has three parts. First, I will build a small media-planning evaluation in R using the vitals package. Then I will discuss why media-planning questions are unusually difficult to evaluate. Finally, I will focus on what the HAI-Q results tell us about general-purpose LLMs.

1. A small eval in R

At its simplest, an LLM eval has three components:

1. A dataset containing questions and targets.
2. A solver that answers the questions.
3. A scorer that evaluates the answers.

This is the basic structure used by vitals. The target column does not always need to contain one exact answer. For open-ended questions, it can contain grading guidance: what a good answer must include, which assumptions must be stated, and which errors should lose credit.

This is useful for testing the final answer of an LLM. It is less useful as a complete test of an AI agent. An agent consists of a harness and one or more models: it may retrieve data, call tools, run calculations, retry a failed step, maintain state, and decide when to ask for help. vitals mainly gives us a dataset, a solver, and a scorer. It does not automatically tell us whether the agent used the right tool, passed the right arguments, or recovered safely from an intermediate failure. It can be extended with custom solvers, but then the agent harness and trajectory checks still have to be built separately.

The questions

Here is a deliberately small evaluation set. The questions cover calculation, planning, and interpretation. The targets are written as instructions for a grader rather than as model answers.

library(dplyr)
library(tibble)

media_questions <- tribble(
  ~id, ~input, ~target, ~domain, ~task,

  "grp-calculation",
  "A campaign delivers 60% reach at an average frequency of 3. What is the campaign's GRP? State the assumption behind the calculation.",
  "A full-credit answer calculates 60 * 3 = 180 GRP. It explains that this is the standard simplified relationship GRP = reach (in percent) * average frequency. It must not describe GRP as the number of unique people reached.",
  "Reach and GRP", "Calculation",

  "budget-allocation",
  "A client has EUR 100,000 to reach adults aged 25-49 in Germany. There is no historical channel performance data. Propose a first-pass media allocation and explain how you would improve it after launch.",
  "There is no single correct allocation. A strong answer makes the uncertainty explicit, explains the role of each channel, avoids inventing CPMs or reach figures, and proposes a measurement and test-and-learn plan. It should distinguish a planning hypothesis from an observed result.",
  "Media strategy", "Recommendation",

  "reach-frequency-tradeoff",
  "Two plans have the same budget. Plan A has higher reach and lower frequency; Plan B has lower reach and higher frequency. How would you decide between them?",
  "A full-credit answer says that the choice depends on the communication objective, the selected channesls, the overall attention, audience size, purchase cycle, creative strength, and expected response. It should explain the reach-frequency trade-off and ask for missing information instead of declaring one plan universally better.",
  "Media strategy", "Reasoning"
)

media_questions

The questions are intentionally specific. “Create a good media plan” is difficult to score because almost any fluent answer can sound reasonable. A question about GRP, in contrast, gives us a concrete numerical check. The budget question tests something different: whether the model knows when it does not have enough information to justify a precise allocation.

In a real eval, I would expand this table with questions about channel roles, target-group definitions, reach curves, media costs, campaign objectives, and German market data. I would also add metadata such as difficulty, source, market, and the date on which the underlying data was valid.

The solver and scorer

The vitals introduction uses Task$new() to combine the dataset, a solver, and a scorer. The following example uses Claude as the solver and an LLM-based grader. An API key for the selected provider is required.

library(vitals)
library(ellmer)

media_task <- Task$new(
  dataset = media_questions,
  solver = generate(
    chat_openrouter(model = "claude-sonnet-4.6")
  ),
  scorer = model_graded_qa(partial_credit = TRUE),
  name = "A small media-planning eval"
)

# Run the solver and scorer.
media_task$eval()

# Bind the results into a tibble for analysis.
media_scores <- vitals_bind(media_task)
media_scores

media_scores |>
  count(score)

With partial_credit = TRUE, the scorer can return C for correct, P for partially correct, or I for incorrect. This is more useful than a single binary score for media planning. A response may use the right calculation but fail to explain its assumptions. Another answer may be directionally useful but invent a media cost. Those are different failures and should not be hidden in one number.

The result is not the end of the evaluation. I would inspect the individual answers and the grading decisions, especially on the first few runs. If the grader repeatedly gives full credit to a response that I consider unsafe, the problem may be the grading guidance rather than the model.

Comparing models

Once the task exists, comparing solvers is straightforward. For example, the same questions can be evaluated with an OpenAI model and then combined with the first run.

media_task_openai <- media_task$clone()

media_task_openai$eval(
  solver_chat = chat_openrouter(model = "gpt-5.4")
)

comparison <- vitals_bind(
  claude = media_task,
  openai = media_task_openai
) |>
  mutate(model = recode(
    task,
    claude = "Claude",
    openai = "OpenAI"
  ))

comparison |>
  count(model, score)

The exact model names will change over time. The important part is to keep the questions and scoring rules stable while changing one component at a time. Otherwise, a score difference is difficult to interpret.

This example is an LLM eval, not an agent eval. It compares the text returned at the end of the solver. For an agent, I would additionally log the tool calls, retrieved documents, intermediate calculations, retries, latency, cost, and final answer. The test should then assert properties of the whole trajectory, not only the final prose.

2. Why media-planning evals are hard

Media planning looks like a natural use case for LLMs. The work involves briefs, audience descriptions, channel recommendations, calculations, and written rationales. But the same combination makes the evaluation difficult.

There is rarely one correct plan

Many planning questions are underdetermined. A sensible allocation depends on the objective, budget, campaign length, audience size, creative assets, buying conditions, historical performance, and measurement setup. Two experienced planners can recommend different channel mixes and both be right.

This means that a target should often describe required properties rather than prescribe one number. For example, a full-credit answer might need to:

• state its assumptions;
• stay within the budget;
• connect channel choices to the objective;
• distinguish facts from hypotheses; and
• explain how the plan would be measured and revised.

The evaluator has to reward sound reasoning without turning one planner’s preference into a false ground truth.

Media facts are local and time-dependent

A media answer can be correct in one market and wrong in another. Reach, costs, inventory, audience definitions, platform availability, and panel methodology all vary by country and by date. A model that has learned general media vocabulary does not automatically have access to the latest German market data.

This is why a useful eval should record the market and the data vintage. It should also test whether the model provides a source or clearly labels an assumption when a current figure is unavailable. “The average CPM is X” is not a safe answer if the model cannot explain where X came from.

Definitions matter more than fluent prose

Media planning has many terms that are close enough to sound interchangeable but are not interchangeable in a calculation. Reach can mean gross or net reach. Frequency can be average or effective frequency. GRP is not the same as the number of unique people reached. Target-group percentages need a denominator.

These distinctions create easy failure modes. A response can contain the right formula but apply it to the wrong quantity. It can also produce a tidy table whose totals do not add up. Numerical checks and constraint checks therefore need to be part of the eval, not just a general assessment of writing quality.

LLM-as-a-judge is useful but not sufficient

The model_graded_qa() scorer in vitals makes open-ended evaluation scalable. It is also a model judging another model’s answer, which introduces its own errors. A grader may prefer confident language, overlook a subtle numerical mistake, or accept an answer because it resembles the supplied guidance.

For a serious benchmark, I would calibrate the grader against human-labelled examples and keep a sample for manual review. I would also separate different dimensions where possible:

• factual accuracy;
• numerical accuracy;
• constraint satisfaction;
• quality of reasoning; and
• usefulness of the recommendation.

One overall score is convenient, but it can hide which part of the planning process is failing.

vitals is not an agent harness

The distinction matters in practice. vitals is a useful way to ask an LLM a fixed set of questions and score the answers. It is not, by itself, a test harness for an AI agent. If an agent has access to a reach database, a calculator, or a planning API, the evaluation needs to check the interaction with those components as well.

For an end-to-end agent eval, I would test at least four layers: whether the harness selects the right tool, whether the tool call is valid, whether the intermediate result is interpreted correctly, and whether the final answer communicates uncertainty. A final response can look reasonable even when one of those earlier steps failed.

3. Results from the Havas HAI-Q

The Havas HAI-Q page describes a benchmark built for media planning in the German market. It contains 35 practical tasks across five areas:

1. Media strategy
2. Channel and budget allocation
3. Target-group evaluation
4. Reach and GRP
5. Media and effectiveness metrics

In this post I focus on the general-purpose LLM results:

The headline is uncomfortable: even GPT-5, the strongest general-purpose model in this comparison, is correct on fewer than half of the questions. It answers 16 of 35 questions correctly, leaving 19 answers that are not counted as correct. For practical media-planning work, that is not a safe error rate.

When I look at the individual HAI-Q questions, the errors are not evenly distributed. The models particularly struggle with numeric questions involving reach figures and audience sizes. They can produce a fluent explanation of a planning idea and still make a basic mistake in a denominator, a percentage, or the interpretation of a reach number. This is exactly the type of failure that is easy to miss if evaluation is based on plausibility rather than calculation.

There is also a clear improvement between the older and newer models. GPT-4o answers 4 of 35 questions correctly, compared with 16 of 35 for GPT-5. That is 12 additional correct answers and an improvement from 11.4% to 45.7%, or 34.3 percentage points. GPT-5 gets four times as many questions correct as GPT-4o on this benchmark. The improvement is substantial, even though the newer model is still not reliable enough for unsupervised media decisions.

Claude Sonnet 4.5 scores lower in this test, with 6 of 35 correct answers. This does not mean that the model is generally poor. It means that a general language benchmark and a domain-specific media-planning benchmark answer different questions. The HAI-Q is designed around the German market and a specific definition of correctness, including the numerical details that matter in planning.

The aggregate scores should still be read with care. They do not show the difficulty of every question, the variance across repeated runs, or the exact scoring decision for every answer. But the pattern is clear enough to be useful: current general-purpose LLMs are good at sounding like media planners, while still struggling to answer roughly half of these practical questions correctly.

Takeaways

An LLM eval for media planning should start with real questions from the workflow, not generic prompts about marketing. In R, a tibble with input and target is enough to get started, and vitals provides a practical way to connect questions, solvers, scorers, and result inspection.

The difficult part is not running the model. It is defining what “correct” means when several plans can be reasonable, data changes by market and date, and a plausible answer can contain a subtle but consequential error. For agents, the evaluation also needs to cover the harness, tools, intermediate steps, and failure handling. A final answer score is not enough.

The HAI-Q results show that the latest general-purpose LLMs still answer more than half of these media-planning questions incorrectly. The largest weaknesses appear in numerical reasoning around reach and audience sizes. At the same time, the jump from GPT-4o to GPT-5 shows that model development is moving in the right direction.

For me, the practical conclusion is simple: use LLMs to accelerate media work, but test every numerical recommendation before relying on it. The next step is to combine model evals such as HAI-Q with end-to-end agent tests that measure tool use and the full planning workflow.

Resources

• Getting started with vitals
• vitals documentation
• HAI-Q: Havas AI Media Quality Index
• ellmer documentation

Like many people I have been keeping track of the slow-burn fuel crisis arising from the war in Iran, its possible expansion to the Red Sea, and exacerbation from the ongoing war between Russia and Ukraine. The latter seems to have reached a stage where serious damage is being done to Russia’s refining capacity and petro products’ exports.

I’ve now started a page on this website with a few charts on the fuel crisis that I’m going to be keeping up to date.

Stocks of crude are going down, but “how much is left”?

The motivation for today’s post was that I’d been seeing many news and opinion articles about how the US crude oil stocks—in total, and in the Strategic Petroleum Reserve in particular—have declined to modern record lows. But expert commentators are reluctant to show this decline in relation to the absolute amount left, as opposed to comparisons to seasonal or historic lows. This is understandable because no-one really knows (or perhaps tells) what the real “tank bottoms” are in a big complex system like the USA. A certain amount of oil, of uncertain volume but certainly in the tens (and possibly hundreds) of millions of barrels, is needed just to keep various bits of plumbing or storage from seizing up or collapsing.

So it’s hard to get a chart like this one, and I had to make it myself:

I’ve got no reputation as an oil expert to worry about, so I can safely draw those lines at a fairly arbitrary 24 (or 30) days of refinery throughput, and calculate the time it will take to get to those low levels “at the current rate”. It’s simplistic, but it’s an important reality check to see that there might be about 19 weeks left.

Fast as the stocks are being used up, the USA has a lot of oil available to it. But 19 weeks isn’t that many. Markets are betting that good things will happen in the weeks or months before the blue line reaches the red line. Because of that, oil is currently selling for around $90 per barrel instead of the $110+ that would be justified if people thought we had lost the traffic through the Strait of Hormuz indefinitely.

Not that I think 19 weeks is a realistic forecast. Long before US stocks reached that red line, there would be economic and political chaos. Efforts to keep retail prices artificially low would probably be overwhelmed and there would be material demand destruction. There would also be huge incentives for the US to give away anything to restore Persian Gulf and Red Sea transit routes. 19 weeks gives us a time frame within which something will have to change.

The value-add for me from this chart is showing the total oil available for the US in the event this becomes a real emergency for them; the comparison of that to a meaningful if arbitrary threshold in terms of days of refinery cover; and the calculation of weeks left at current draw-down rates to reach that high stress threshold.

Break out the SPR and also look at gasoline, distillates

The same data source has a lot of other detailed information, most of it far too detailed for me, but I’m interested enough in four key headline variables that are more detailed than “total stocks of crude oil”. In this chart I break down the total stocks into those in the Strategic Petroleum Reserve, and those in more normally accessed stocks. I also show the amounts held of two key refined products:

Code

Here’s the code that grabs the data from the US Energy Information Administration and draws that second plot:

library(tidyverse)
library(readxl)
library(scales)
library(glue)

fn <- "us-petrol-status-weekly.xls"

download.file(
  "https://ir.eia.gov/wpsr/psw01.xls",
  fn,
  mode = "wb"
)

us_stocks <- read_excel(fn, sheet = "Data 1", skip = 2) |>
  rename(
    crude = `Weekly U.S. Ending Stocks of Crude Oil  (Thousand Barrels)`,
    crude_spr = `Weekly U.S. Ending Stocks of Crude Oil in SPR  (Thousand Barrels)`,
    gasoline = `Weekly U.S. Ending Stocks of Total Gasoline  (Thousand Barrels)`,
    diesel = `Weekly U.S. Ending Stocks of Distillate Fuel Oil  (Thousand Barrels)`,
    date = Date
  )

#------------------plotting------------------

eia_caption <- glue("Source: Energy Information Administration (EIA). Accessed {format(Sys.Date(), '%d %B %Y')}")

# an annotation rectangle we can use in multiple charts
war_rect <- annotate(
  "rect",
  xmin = as.Date("2026-02-28"),
  xmax = Inf,
  ymin = -Inf,
  ymax = Inf,
  alpha = 0.5,
  fill = "grey80"
)

# Facet plot of total crude, SPR crude, diesel and gasoline stocks
us_stocks |>
  mutate(com_crude = crude - crude_spr) |>
  select(
    date,
    `Non-SPR crude oil` = com_crude,
    `Crude oil in SPR` = crude_spr,
    `Gasoline` = gasoline,
    `Distillate (mostly diesel)` = diesel
  ) |>
  gather(variable, value, -date) |>
  mutate(variable = fct_reorder(variable, value)) |>
  ggplot(aes(x = date, y = value / 1000)) +
  war_rect +
  facet_wrap(~variable, scales = "free_y") +
  geom_line(colour = "blue") +
  expand_limits(y = 0) +
  scale_y_continuous(label = comma) +
  labs(
    x = "",
    y = "Thousands of barrels",
    title = "US stocks of crude oil, gasoline, and distillate fuel oil (effectively diesel)",
    subtitle = "Showing both total crude oil stocks (crude) and those in the Strategic Petroleum Reserve (crude_spr)",
    caption = eia_caption
  )

As I side point, I’ve started using air to format my code. Nearly always it comes out looking better and easier to read than my manual formatting, and it’s certainly more consistent. Thanks to D. Vaughan and L. Henry and the tidyverse.org project.

And here’s the code to draw my highly polished main plot, the one of total crude oil stocks with the annotations. I’ve tried to do this in a way that it will keep working in the months ahead with minimal maintenance, we’ll see how that goes.

refinery_throughput <- 17.3 # as at 24 July, operating at 97% of US capacity. Daily refinery use.
plausible_stress <- 30 * refinery_throughput
plausible_high_stress <- 24 * refinery_throughput

crude_growth_summary <- us_stocks |>
  summarise(
    latest_crude = crude[date == max(date)],
    weeks = as.numeric(as.Date(max(date)) - as.Date("2026-04-03")) / 7,
    difference = (latest_crude -
      crude[date == as.Date("2026-04-03")]) /
      1000,
    ratio = latest_crude / crude[date == as.Date("2026-04-03")],
    # growth/decline per week in million barrels e.g. 10m barrels per week:
    difference_rate = difference / weeks,
    # growth/declien rate per week:
    growth_rate = 1 - exp(log(ratio) / weeks)
  ) |>
  mutate(
    weeks_at_this_rate = (plausible_stress - latest_crude / 1000) /
      difference_rate
  )

plot_us_stocks <- function(min_date = "2020-01-01", lab_x_diff = NULL) {
  if (!"Date" %in% class(min_date)) {
    min_date <- as.Date(min_date)
  }

  if (is.null(lab_x_diff)) {
    days_shown <- as.numeric(as.Date(max(us_stocks$date))) -
      as.numeric(min_date)
    lab_x_diff <- days_shown / 50
  }

  lv <- tail(us_stocks, 1)$crude / 1000

  p2 <- us_stocks |>
    filter(date >= min_date) |>
    ggplot(aes(x = date, y = crude / 1000)) +
    war_rect +
    geom_hline(yintercept = plausible_stress, colour = "darkred") +
    geom_hline(yintercept = plausible_high_stress, colour = "red") +
    annotate(
      "text",
      x = min_date + lab_x_diff,
      y = plausible_stress + 50,
      label = "Illustrative stress threshold:\n30 days of refinery cover",
      colour = "darkred",
      hjust = 0
    ) +
    annotate(
      "text",
      x = min_date + lab_x_diff,
      y = plausible_high_stress - 30,
      label = "Illustrative high stress threshold:\n24 days of refinery cover",
      colour = "red",
      hjust = 0,
      vjust = 1
    ) +
    annotate(
      "text",
      x = max(us_stocks$date) - 5e6,
      y = lv,
      label = glue(
        "{round(lv)} million barrels;\n{round(lv / refinery_throughput)} days of cover"
      ),
      size = 2.9,
      hjust = 1,
      vjust = 1,
      colour = "blue"
    ) +
    geom_line(colour = "blue") +
    expand_limits(y = 0) +
    scale_y_continuous(
      label = comma,
      sec.axis = sec_axis(
        ~ . / refinery_throughput,
        name = "Days of refinery throughput"
      )
    ) +
    labs(
      x = "",
      y = "Millions of barrels",
      title = "U.S. Total Crude Oil Stocks, including Strategic Petroleum Reserve",
      subtitle = glue(
        "Comparison of existing inventories with refinery throughput as at July 2026.
Decline since peak on 3 April 2026 is at {abs(round(crude_growth_summary$difference_rate, 1))} million barrels per week; {round(crude_growth_summary$weeks_at_this_rate)} weeks from stress threshold at this (hypothetical and linear) rate."
      ),
      caption = eia_caption
    ) +
    theme(axis.line.y = element_line(colour = "grey50"))

    svglite(
      glue("all-us-crude-from{min_date}.svg"),
      width = 10,
      height = 7
    )
    print(p2)
    dev.off()
}

plot_us_stocks("2000-01-01")
plot_us_stocks("2020-01-01")
plot_us_stocks("2025-01-01")

Because I did this plot in the form of a function, I have a couple of variants of the plot. One starting in 2020:

And the other zoomed right in to just January 2025 and onwards:

On the whole I prefer the plot I started with, showing the situation from 2000 and onwards. I’m trying to get some big picture perspective after all.

That’s all for today. The images on this page will stay the same as at the time of writing; but those on the fuel crisis monitoring page will be updated each week.

era v0.6.0 is now available on CRAN:

install.packages("era")
library(era)

This minor release adds functions for chronological comparison of year vectors:

x <- yr(c(200, 100, 300), "BCE")
yr_earlier_than(x, yr(150, "BCE"))
#> # BCE years <yr_lgl[3]>:
#> [1] FALSE FALSE  TRUE
#> # Era: Before Common Era (BCE): Gregorian years (365.2425 days), counted backwards from 1

yr_later_than(x, yr(150, "BCE"))
#> # BCE years <yr_lgl[3]>:
#> [1]  TRUE FALSE FALSE
#> # Era: Before Common Era (BCE): Gregorian years (365.2425 days), counted backwards from 1

And for calculating chronological differences between year vectors:

yr_difference(yr(300, "BCE"), yr(100, "BCE"))
#> # years <yr[1]>:
#> [1] 200
#> # Era: Gregorian years (365.2425 days)

These functions are aware of era directionality: “BCE” years count backwards, so 300 BCE is earlier than 100 BCE. They support work on tempo, a new package I am working on that represents chronological intervals and the temporal relations between them.

Links

• era on CRAN
• era package documentation
• Introductory vignette
• Source code (GitHub)

The goal of spuriouscorrelations is to keep alive the amazing examples from Tyler Vigen. Unfortunately, as of 2023-10-09, the website is down as my students noticed. Therefore, I decided to use the snapshot from the Internet Wayback Machine to save the datasets from 2023-06-07.

Installation

You can install the CRAN version of spuriouscorrelations with:

install.packages("spuriouscorrelations")

You can install the development version of spuriouscorrelations with:

remotes::install_github("pachadotdev/spuriouscorrelations")

Example

The package covers multiple examples for different spurious (and curious) correlations.

Here is a basic example which shows you how to plot a spurious correlation for the variables

• x: Number of people who drowned by falling into a pool
• y: Films Nicolas Cage appeared in

library(spuriouscorrelations)
library(tinyplot)

pool_drownings
   year   x y
1  1999 109 2
2  2000 102 2
3  2001 102 2
4  2002  98 3
5  2003  85 1
6  2004  95 1
7  2005  96 2
8  2006  98 3
9  2007 123 4
10 2008  94 1
11 2009 102 4
cor(pool_drownings$x, pool_drownings$y)
[1] 0.6660043
tinyplot(
  y ~ x,
  data = pool_drownings,
  main = sprintf("Correlation %s", round(cor(pool_drownings$x, pool_drownings$y), 3)),
  xlab = "Number of people who drowned by falling into a pool",
  ylab = "Films Nicolas Cage appeared in"
)

Converting the data to long format simplifies plotting both variables per year:

pool_drownings_2 <- reshape(
  pool_drownings, 
  varying = c("x", "y"),  # columns to collapse
  v.names = "value",      # name of the new value column
  timevar = "variable",   # name of the new ID column
  times = c("x", "y"),    # values to populate the ID column
  direction = "long"      # target format
)

pool_drownings_2
     year variable value id
1.x  1999        x   109  1
2.x  2000        x   102  2
3.x  2001        x   102  3
4.x  2002        x    98  4
5.x  2003        x    85  5
6.x  2004        x    95  6
7.x  2005        x    96  7
8.x  2006        x    98  8
9.x  2007        x   123  9
10.x 2008        x    94 10
11.x 2009        x   102 11
1.y  1999        y     2  1
2.y  2000        y     2  2
3.y  2001        y     2  3
4.y  2002        y     3  4
5.y  2003        y     1  5
6.y  2004        y     1  6
7.y  2005        y     2  7
8.y  2006        y     3  8
9.y  2007        y     4  9
10.y 2008        y     1 10
11.y 2009        y     4 11
tinyplot(
  value ~ year | variable, # "|" indicates the groping variable for the legend
  data = pool_drownings_2,
  main = sprintf("Correlation %s", round(cor(pool_drownings$x, pool_drownings$y), 3)),
  xlab = "Year",
  ylab = "Pooled observations",
  pch = 19 # solid dot shape
)

How about standarzing the variables to avoid the different scale visibility issue?

pool_drownings_3 <- pool_drownings

pool_drownings_3$x <- (pool_drownings_3$x - mean(pool_drownings_3$x)) / sd(pool_drownings_3$x)
pool_drownings_3$y <- (pool_drownings_3$y - mean(pool_drownings_3$y)) / sd(pool_drownings_3$y)

pool_drownings_3 <- reshape(
  pool_drownings_3, 
  varying = c("x", "y"),  # columns to collapse
  v.names = "value",      # name of the new value column
  timevar = "variable",   # name of the new ID column
  times = c("x", "y"),    # values to populate the ID column
  direction = "long"      # target format
)

pool_drownings_3
     year variable      value id
1.x  1999        x  0.8952867  1
2.x  2000        x  0.1696333  2
3.x  2001        x  0.1696333  3
4.x  2002        x -0.2450258  4
5.x  2003        x -1.5926679  5
6.x  2004        x -0.5560201  6
7.x  2005        x -0.4523554  7
8.x  2006        x -0.2450258  8
9.x  2007        x  2.3465935  9
10.x 2008        x -0.6596849 10
11.x 2009        x  0.1696333 11
1.y  1999        y -0.2470999  1
2.y  2000        y -0.2470999  2
3.y  2001        y -0.2470999  3
4.y  2002        y  0.6589330  4
5.y  2003        y -1.1531327  5
6.y  2004        y -1.1531327  6
7.y  2005        y -0.2470999  7
8.y  2006        y  0.6589330  8
9.y  2007        y  1.5649658  9
10.y 2008        y -1.1531327 10
11.y 2009        y  1.5649658 11
tinyplot(
  value ~ year | variable, # "|" indicates the groping variable for the legend
  data = pool_drownings_3,
  main = sprintf("Correlation %s", round(cor(pool_drownings$x, pool_drownings$y), 3)),
  xlab = "Year",
  ylab = "Pooled standardized observations",
  pch = 19, # solid dot shape
  type = "b" # add line to see the trend clearly
)

Join our workshop on Introduction to All Things About bpvars, the R package for Forecasting with Bayesian Panel Vector Autoregressions,  which is a part of our workshops for Ukraine series! 

Here’s some more info: 

Title: All Things About bpvars, the R package for Forecasting with Bayesian Panel Vector Autoregressions

Date: Thursday, September 10th, 18:00 – 20:00 CEST (Rome, Berlin, Paris timezone) 

Speaker: Tomasz Woźniak is a Bayesian econometrician developing new econometric methods for applied macroeconomic research. He has been a specialised R user for eighteen years, and has recently joined The R Journal as an Associate Editor. He is the author of several R packages, available at https://bsvars.org/, that combine blazingly fast algorithms written in C++ with the convenience of data analysis in R. He works as a senior lecturer at the University of Melbourne, where he has an extensive research, teaching, and engagement portfolio. Tomasz supports Ukraine. Slava Ukrainie!

Description: This session is a unique opportunity to learn to use the bpvars package for forecasting with Bayesian panel vector autoregressions in a wide range of applications spanning:

* data preparation,

* specification of panel vector autoregressions,

* estimation of the models,

* estimation with missing observations,

* forecasting labour market outcomes globally,

* plotting and reporting predictions of different qualities, including point and density forecasts, marginal, conditional, or restricted forecasts,

* performing recursive expanding window forecasting,

* generating forecasting performance reports.

A sequence of hands-on exercises that fully prepares attendees to work with the package supports all of this. You are welcome to browse existing resources, including:

* bpvars website: https://bsvars.org/bpvars/

* bpvars CRAN profile: https://cran.r-project.org/package=bpvars

* package vignette: https://doi.org/10.48550/arXiv.2606.14143

* other materials: https://bsvars.org/bpvars/#resources

See you soon!

Preparation:

This session assumes the attendees know the basics of time series analysis. Install and load the bpvars package, then check the documentation by running the code below. Run the example you will find there in R. If it works, you’re ready to go. Thanks!

“`

install.packages(“bpvars”)

library(bpvars)

?bpvars

Minimal registration fee: 20 euro (or 20 USD or 800 UAH)

Please note that the registration confirmation is sent 1 day before the workshop to all registered participants rather than immediately after registration

How can I register?

• Go to https://bit.ly/3wvwMA6 or https://bit.ly/4aD5LMC  or  https://bit.ly/3PFxtNA and donate at least 20 euro. Feel free to donate more if you can, all proceeds go directly to support Ukraine.

• Save your donation receipt (after the donation is processed, there is an option to enter your email address on the website to which the donation receipt is sent)

• Fill in the registration form, attaching a screenshot of a donation receipt (please attach the screenshot of the donation receipt that was emailed to you rather than the page you see after donation).

If you are not personally interested in attending, you can also contribute by sponsoring a participation of a student, who will then be able to participate for free. If you choose to sponsor a student, all proceeds will also go directly to organisations working in Ukraine. You can either sponsor a particular student or you can leave it up to us so that we can allocate the sponsored place to students who have signed up for the waiting list.

How can I sponsor a student?

• Go to https://bit.ly/3wvwMA6 or https://bit.ly/4aD5LMC  or https://bit.ly/3PFxtNA and donate at least 20 euro (or 17 GBP or 20 USD or 800 UAH). Feel free to donate more if you can, all proceeds go to support Ukraine!

• Save your donation receipt (after the donation is processed, there is an option to enter your email address on the website to which the donation receipt is sent)

• Fill in the sponsorship form, attaching the screenshot of the donation receipt (please attach the screenshot of the donation receipt that was emailed to you rather than the page you see after the donation). You can indicate whether you want to sponsor a particular student or we can allocate this spot ourselves to the students from the waiting list. You can also indicate whether you prefer us to prioritize students from developing countries when assigning place(s) that you sponsored.

If you are a university student and cannot afford the registration fee, you can also sign up for the waiting list here. (Note that you are not guaranteed to participate by signing up for the waiting list).

You can also find more information about this workshop series,  a schedule of our future workshops as well as a list of our past workshops which you can get the recordings & materials here.

Looking forward to seeing you during the workshop!

I’m a big fan of Github Copilots ghost text and next edit suggestions. As a data analyst they give you much more control over the pace and direction of analysis code than a full AI agent. They also let you practice the core skill an agent needs from you: writing a clear specification.

I’ll walk through it in R, but the same ideas apply in any language. I’m assuming you’re on VScode with the Copilot extension installed. This should work with the free or paid plans (though I haven’t checked free myself).

library(dplyr)
library(ggplot2)
dat <- data.frame(
    time = rep(1:10, 3),
    species = rep(c("A", "B", "C"), each = 10),
    value = c(rnorm(10, mean = 5), rnorm(10, mean = 10), rnorm(10, mean = 15))
)

datA <- dat |> filter(species == "A")
m1 <- lm(value ~ time, data = datA)
coef(m1)["time"]

ggplot(datA, aes(x = time, y = value)) +
    geom_point() +
    geom_smooth(method = "lm", se = FALSE) +
    labs(
        title = "Species A Time Series",
        x = "Time",
        y = "Value"
    )

fit_fun <- function(species_name) {
    datA <- dat |> filter(species == species_name)
    m1 <- lm(value ~ time, data = datA)
    coef(m1)["time"]

    ggplot(datA, aes(x = time, y = value)) +
        geom_point() +
        geom_smooth(method = "lm", se = FALSE) +
        labs(
            title = paste("Species", species_name, "Time Series"),
            x = "Time",
            y = "Value"
        )
}

fit_fun("A")

# Simulate a new dataset of abundance at x-y coordinates

# plot a 2D map

# fit a model with interaction between x and y

Four hundred twelve new packages were submitted to CRAN in June. Here are my Top 40 picks in nineteen categories: Bioarchaeology, Biology, Climate Studies, Computational Methods, Ecology, Epidemiology, Finance, Functional Data Analysis, Machine Learning, Medical Statistics, Networks, Pharmacokinetics, Probability, Programming, Psychometrics, Risk Analysis, Statistics, Time Series, and Utilities.

This post describes the R version of Python’s GPopt (https://docs.techtonique.net/GPopt/GPopt.html), a package for ‘Bayesian’ optimization of black-box functions and machine learning hyperparameter tuning, using Gaussian Process Regression and other conformalized surrogates. The package is available on GitHub, and through the R universe.

Keep in mind that this package is for Machine Learning hyperparameter tuning: the global minimum won’t always be found, but this isn’t an issue, since it means you aren’t overfitting the training set.

It’s ported the same way as nnetsauce for R was: with uv to create an isolated Python virtual environment containing the Python GPopt package, and reticulate to call into it from R. Every function in this R package is a thin wrapper that returns the underlying Python object; the general rule is: object accesses with .’s in Python are replaced by $’s in R.

See this post for the technique: Finally figured out a way to port python packages to R using uv and reticulate.

Install

1. Create a Python virtual environment with uv

# pip install uv # if necessary
uv venv venv
source venv/bin/activate      # on Windows: venv\Scripts\activate
uv pip install pip GPopt

Keep track of where venv/ lives – you’ll pass its path as venv_path to every function in this package.

2. Install the R package

install.packages("remotes")
remotes::install_github("Techtonique/GPopt_r") 

reticulate will be installed automatically as a dependency.

Examples

Minimizing the Branin function

This is a standard test function for optimization algorithms. GPOpt is more suitable for expensive black-box functions, but this is a good example to illustrate the usage of the package.

library(GPopt)

branin <- function(x) {
  x1 <- x[1]; x2 <- x[2]
  term1 <- (x2 - (5.1 * x1^2) / (4 * pi^2) + (5 * x1) / pi - 6)^2
  term2 <- 10 * (1 - 1 / (8 * pi)) * cos(x1)
  term1 + term2 + 10
}

opt <- GPOpt(
  lower_bound = c(-5, 0),
  upper_bound = c(10, 15),
  objective_func = branin,
  n_init = 10,
  n_iter = 40,
  venv_path = "./venv"
)

opt$optimize(verbose = 1L)
print(opt$x_min)  # best parameters
print(opt$y_min)  # best objective value

Tuning a scikit-learn model’s hyperparameters

library(GPopt)

sklearn <- get_sklearn(venv_path = "./venv")
RandomForestClassifier <- sklearn$ensemble$RandomForestClassifier

X <- as.matrix(iris[, 1:4])
y <- as.integer(iris$Species) - 1L

mlopt <- MLOptimizer(scoring = "accuracy", cv = 5, venv_path = "./venv")

param_config <- list(
  n_estimators = list(bounds = c(10, 300), dtype = "int"),
  max_depth    = list(bounds = c(1, 20),   dtype = "int")
)

mlopt$optimize(
  X_train = X, y_train = y,
  estimator_class = RandomForestClassifier(),
  param_config = param_config,
  verbose = 1L
)

print(mlopt$get_best_parameters())
print(mlopt$get_best_score())

Bayesian optimization with early stopping

library(GPopt)

opt <- BOstopping(
  f = branin,
  bounds = rbind(c(-5, 10), c(0, 15)),
  venv_path = "./venv"
)
result <- opt$optimize(n_iter = 100L)

Using a custom (conformalized) surrogate model

library(GPopt)

sklearn <- get_sklearn(venv_path = "./venv")
ns <- get_nnetsauce(venv_path = "./venv")

opt <- GPOpt(
  lower_bound = c(-5, 0),
  upper_bound = c(10, 15),
  objective_func = branin,
  acquisition="ucb",
  method="splitconformal",
  surrogate_obj = ns$PredictionInterval(sklearn$ensemble$RandomForestRegressor()),
  venv_path = "./venv"
)
opt$optimize(verbose = 1L)

Biological phenomena are often studied by examining how a numerical variable, usually called the response (e.g., the weight or height of an organism), is affected by another variable, usually called the predictor (e.g., time or fertiliser application rate). These ‘predictor-response’ relationships are commonly described using models, which express the response as a mathematical function of the predictor. In general, a model can be written as:

where is the response, is the predictor, and is a collection of parameters, usually denoted by symbols, such as the letters of the Roman or Greek alphabet. The element is the function, which determines the shape of the relationship when plotted on an x–y graph.

Because mathematical modelling plays such a central role in biology and many other scientific disciplines, biologists need to be familiar with the most important mathematical functions. More importantly, they need to be able to “read” these functions and use their parameters to describe, interpret, and quantify biological processes. With this aim in mind, I have compiled a collection of the mathematical functions most commonly encountered in biology, explaining the meaning of their parameters, with particular emphasis on their biological interpretation rather than their mathematical properties.

Curve shapes

Functions are often classified according to the shape they exhibit when plotted on an x–y graph. This approach is adopted, for example, in Ratkowsky (1990), and the classification presented below is largely based on that work.

• Polynomials
  1. Straight line function
  2. Quadratic polynomial function
• Concave/Convex curves (no inflection)
  1. Exponential function
  2. Asymptotic function / Negative exponential function
  3. Power function
  4. Logarithmic function
  5. Rectangular hyperbola
• Sigmoidal curves
  1. Logistic function
  2. Gompertz function
  3. Modified Gompertz function
  4. Log-logistic function
  5. Weibull (type-1) function
  6. Weibull (type-2) function
• Curves with maxima/minima
  1. Peaked sigmoidal function
  2. Bragg function
  3. Lorentz function
  4. Beta function

I have chosen to include a relatively large number of functions, so this post is necessarily rather long. To make navigation easier, you can first inspect Figure 1 to identify the type of curve of interest and then use the links above to jump directly to the corresponding section.

Throughout this post, denotes the predictor, denotes the response, and the letters of the Roman alphabet (, , , …) denote the model parameters, which determine the shape of the function.

D(expr = expression(a + b*X + c*X^2), name = "X")

b + c * (2 * X)

Concave/Convex curves (no inflection points)

Exponential function

The most common parameterisation of an exponential function is:

Depending on the sign of , Eq. 3 describes a concave-up monotonically increasing shape (exponential growth; ) or a concave-up monotonically decreasing shape (exponential decay; ).

The parameter represents the value of when , whereas the meaning of can be understood by considering the first derivative:

D(expression(a * exp(k * X)), "X")

a * (exp(k * X) * k)

In algebraic terms, it is:

and thus:

Therefore, represents the relative rate of change (increase/decrease), throughout the domain, which is often known in growth analysis as the Relative Growth Rate (RGR).

Other equivalent parameterisations are:

and:

Equations 3-5 are equivalent, as can be shown by setting and :

Another, slightly different, parameterisation is commonly used in bioassay studies, mainly to describe exponential decay processes:

where corresponds to in the previous models and .

For all the exponential decay equations presented above, as . A lower asymptote, , can be incorporated into Eq. 6 to account for situations in which the response does not approach zero as the predictor tends to infinity:

Exponential functions are often used to describe the growth of populations under non-limiting environmental conditions or the degradation of xenobiotics in the environment (first-order degradation kinetics). In both cases, represents time and is therefore restricted to non-negative values, while (or, equivalently, ) must be different from zero. Check the post at this link for an example of fitting.

Asymptotic function

This function appears in several different parameterisations and is also known as the monomolecular growth function, the Mitscherlich law, or the von Bertalanffy law. Owing to its biological interpretation, the most widely used parameterisation is:

It describes a monotonically increasing, concave-down curve approaching a horizontal asymptote as tends to infinity. The parameter represents the maximum attainable value of (the plateau), while is the value of when (the initial value). The parameter is proportional to the rate at which the response approaches the plateau. Indeed, the first derivative is

D(expression(a - (a - b) * exp (- m * X)), "X")

(a - b) * (exp(-m * X) * m)

Considering Eq.8, we can write:

and then:

The previous expression shows that the relative rate of change in the response (the RGR in growth analysis) is not constant, as it is for the exponential function. Instead, it depends on the attained value of . In particular, the RGR is greatest at the beginning of the process, when is smallest, and gradually approaches zero as approaches the plateau .

Another closely related parameterisation is often encountered, in which is replaced by , by , and by . Some simple algebraic manipulations are also introduced to make the model consistent with the parameterisation commonly adopted for biological assay models:

In some applications, is reparameterised through its logarithm to facilitate model fitting:

For all these equations, setting (or equivalently ) yields a curve passing through the origin. In this case, the most common parameterisation is:

which is known as the negative exponential function.

The asymptotic function is widely used in biology, not only for growth analysis but also, for example, in weed competition studies. The negative exponential function has also been used to model the absorbed Photosynthetically Active Radiation () as a function of leaf area index (). In this context, represents the incident PAR (), while is the light extinction coefficient. Check the post at this link for an example of fitting.

Power function

The power function is also known as Freundlich function or allometric function. The most common parameterisation is:

Its shape is highly flexible and is determined by the value of the parameter (see Figure 1). When , the response increases with and the curve is concave down. When , decreases as increases and the curve is concave up. Finally, when , increases with and the curve is concave up. These three cases are illustrated in Figure 1 using different colours. The function is defined only for and has no horizontal asymptote as .

The biological interpretation of the parameters and is not straightforward. Both influence the slope of the curve, as can be seen by examining the first derivative:

D(expression(a * X^b), "X")

a * (X^(b - 1) * b)

The power function is mathematically equivalent to an exponential function of , since:

The power function (named as the Freundlich equation) is widely used in agricultural chemistry, for example to model the sorption of xenobiotics in soil. It is also used to describe the relationship between the number of plant species and the sampled area (the species–area relationship). Check the post at this link for an example of fitting.

Logarithmic function

This function is linear in :

Because of the logarithmic transformation, the function is only defined for . The parameter determines the shape of the curve. When , the response increases with and the curve is concave down. Conversely, when , decreases as increases and the curve is concave up, as shown in Figure 1.

The interpretation of the parameters is fairly straightforward. The parameter is the value of the response when . Changing while keeping constant simply shifts the curve vertically without altering its shape, so that the difference between any two curves remains constant for every value of (Fig. 2).

The parameter controls the slope of the curve and, more specifically, it is equal to the slope at , as can be seen from the expression for the first derivative:

D(expression(a + b*log(X)), "X")

b * (1/X)

Changing while keeping constant produces curves that all intersect at (Fig. 2).

In biology, logarithmic functions are used, for example, to describe species–area relationships in ecology, enzyme kinetics in biochemistry, and sensory perception in neurobiology. Check the post at this link for an example of fitting.

Rectangular hyperbola

The rectangular hyperbola is commonly known as the Michaelis-Menten function and is usually parameterised as

The curve is monotonically increasing and concave down, approaching the horizontal asymptote (plateau) . It passes through the origin of the axes ( when ). The parameter is the value of that produces a response equal to . Indeed,

which is readily solved to obtain .

The first derivative is:

D(expression( (a*X) / (b + X) ), "X")

a/(b + X) - (a * X)/(b + X)^2

From this expression, we can see that the initial slope (at ) is .

Equation 14 is not defined for and has no biological meaning for . In most biological applications, however, both and are positive.

An equivalent parameterisation is:

In this form, the parameters and are often replaced by and , respectively. Although the response is not defined for , it tends to zero as .

Another common parameterisation includes the initial slope as an explicit parameter because of its biological relevance. This is obtained by dividing both the numerator and denominator by and noting that , so that :

The Michaelis-Menten function is widely used in pesticide chemistry, enzyme kinetics (Eq. 14), biological assay models (Eq. 15), and weed competition studies, where it is used to describe crop yield loss as a function of weed density (Eq. 16). An example of fitting is shown in another post, at this link.

Sigmoidal functions

Sigmoidal functions are S-shaped, with two horizontal asymptotes (a lower and an upper asymptote) and an inflection point. They can be parameterised in countless ways, which may be confusing when selecting the most appropriate model for a particular biological process. In this post, I will use the parameterisations proposed in Ritz et al (2019), which is very consistent and puts all the curves on the same ground. Furthermore, the fitting properties are good, which contribute to good convergence and reliable estimation.

For the sake of simplicity, it is useful to recognise that most sigmoidal models belong to one of three basic families, which differ only in the position of the inflection point relative to the two asymptotes:

1. Logistic: the ordinate of the inflection point lies halfway between the lower and upper asymptotes; consequently, the S-shape is symmetric (the concave-up and the concave-down portions ‘mirror’ each other).
2. Gompertz: the ordinate of the inflection point lies closer to the lower asymptote, and, thus, the concave-up portions is ‘shorter’ than the concave-down portion.
3. Modified Gompertz: the ordinate of the inflection point lies closer to the upper asymptote, and, thus, the concave-up portion is ‘longer’ than than the concave-down portion.

Each of these three equations can describe either increasing or decreasing responses, depending on the sign of one of the model parameters (see below). Their different shapes are compared side by side in Fig. 3.

Each of these three families (logistic, Gompertz and modified Gompertz) can be transformed by replacing with , giving raise to three corresponding families:

1. Log-logistic
2. Type-I Weibull
3. Type-II Weibull

Like their linear-scale counterparts, these three families can describe either increasing or decreasing responses, depending on the sign of the same model parameter (see below). In the end, we have twelve different sigmoidal functions to choose from. The choice among them depends primarily on the expected symmetry of the response and on whether the predictor is more naturally interpreted on a linear or logarithmic scale.

Thanks to their versatility, sigmoidal functions have been used for an uncountable number of biological applications. Just to mention a few examples in agriculture, I would like to cite plant growth, dose-response curves in biological assays and the time-course of seed germination. For some examples, take a look at my post at this link.

Logistic function

The logistic curve is derived from the cumulative logistic distribution function. It is symmetric about its inflection point and can be parameterised as:

where is the upper asymptote, is the lower asymptote, is the value of at the inflection point, and is proportional to the slope at the inflection point. Because the curve is symmetric, also represents the value of that produces a response halfway between and (commonly referred to as the ED50 in biological assays). The parameter may be either positive or negative and, consequently, the response may either increase or decrease as increases.

This equation is known as the four-parameter logistic model. If appropriate, constraints can be imposed on the parameter values. For example, may be fixed at 0, yielding the three-parameter logistic model. If, in addition, is fixed at 1, the model reduces to the two-parameter logistic model.

Gompertz function

The Gompertz curve can be parameterised in many different ways. I prefer a parameterisation that closely resembles that of the logistic function:

Unlike the logistic function, the Gompertz curve is not symmetric about its inflection point. It exhibits a longer lag phase at the beginning, followed by a progressively steeper increase before gradually approaching the upper asymptote. The parameters have essentially the same interpretation as those of the logistic function, except that , the abscissa of the inflection point, does not correspond to the value of producing a response halfway between and .

As with the logistic function, four-, three-, and two-parameter Gompertz models can be obtained by constraining one or both asymptotes.

Modified Gompertz function

We have seen that the logistic curve is symmetric about its inflection point, whereas the Gompertz curve exhibits a longer lag phase at the beginning, followed by a progressively steeper increase. A different asymmetric pattern can be obtained by modifying the Gompertz function as follows:

The resulting curve increases rapidly at the beginning but gradually slows down as it approaches the upper asymptote. As with the logistic and Gompertz functions, one or both asymptotes can be constrained ( and/or ), giving rise to four-, three-, and two-parameter versions of the modified Gompertz model.

Log-logistic function

The log-logistic curve is symmetric with respect to . A log-normal curve has a very similar shape, although it is used much less frequently. In biological assays (and also in germination studies), the log-logistic function is commonly parameterised as

The parameters have the same interpretation as in the logistic function. In particular, represents the value of that produces a response halfway between and (the ED50). It is easy to show that the above equation is equivalent to

Like the logistic function, the log-logistic model can be fitted in four-, three-, or two-parameter versions by constraining one or both asymptotes. It is widely used to describe dose-response relationships in biological assays, seed germination, and crop growth.

Weibull function (type 1)

The Type I Weibull function is the logarithmic counterpart of the Gompertz function, being defined on rather than on . It is parameterised as

The parameters have essentially the same interpretation as those of the other sigmoidal functions presented above. In particular, and are the lower and upper asymptotes, respectively, while is the value of corresponding to the inflection point. Unlike the log-logistic function, however, does not correspond to the ED50.

Weibull function (type 2)

The Type II Weibull function is closely related to the Type I Weibull function, but it describes a different type of asymmetry, analogous to that of the modified Gompertz function:

The parameters have the same interpretation as those of the Type I Weibull function. In particular, and are the lower and upper asymptotes, respectively, while is the value of at the inflection point. As with the Type I Weibull function, does not correspond to the ED50. One or both asymptotes may be constrained, giving rise to four-, three-, and two-parameter versions of the model.

Another ‘flexible’ sigmoid

I would also like to mention another sigmoidal function that has been widely used in biology because of its supposed flexibility: the Richards function:

In this model, the degree of asymmetry is controlled by the parameter . When , the function reduces to the logistic model, whereas, as decreases below 1, the curve progressively resembles a Gompertz function. Figure 4 illustrates the different shapes obtained for different values of .

Despite this apparent flexibility, the Richards function has poor statistical properties for parameter estimation. In particular, the additional parameter often causes strong correlations among the estimates, making the fitting process unstable and the parameter estimates difficult to interpret. Ratkowsky (1990) described it as having “more undesirable nonlinear regression behaviour than almost any nonlinear regression model in common use.” For these reasons, the Richards function will not be considered further in this blog.

      As a kid, my baseball hero was Mickey Mantle. In his 18-year Major League career, he amassed tremendous statistics, and surely he was one of the game’s all-time greats. When Mantle retired, he was third on the all-time career home run list with 536 home runs, trailing only Babe Ruth (714) and Willie Mays (587). Unfortunately, he suffered numerous injuries, leaving fans to wonder what his numbers would look like if he had enjoyed a healthier career.

      For decades many people considered (and still consider) Babe Ruth to be the greatest all-around baseball player, because he was a star pitcher and a star hitter. The counter-argument to the Mantle “what-if” is: what if Ruth had not spent his first six years as a pitcher and instead had been a full-time hitter? (Similarly, what if Ted Williams had not spent some of his prime years in the military.)

      Nowadays Shohei Ohtani is certainly a once-in-a-generation two-way star. However, he spent his early prime years in Japan’s professional baseball, and he may not have enough time left in US baseball to amass these historic lifetime stats. Aaron Judge is another exceptional modern hitter; hopefully, he stays healthy so we can watch his career numbers climb.

      To compare Mantle and Ruth’s lifetime stats objectively, I turned to data. The Sean Lahman baseball dataset contains Major League player stats back to 1871, and is available in the R library lahman. This was a good opportunity to practice data manipulation using dplyr.

      Ruth played 102 more games than Mantle over a 22-year career (2,503 versus 2,401). His lifetime statistics eclipse Mantle’s in every offensive category except stolen bases — though Ruth did steal home 10 times and hit nearly twice as many triples, suggesting he was faster than most realize. His lifetime batting average was an impressive .342, sitting just behind Tris Speaker and Ted Williams.

      Lineup protection played a massive role for both men. Ruth generally batted directly before Lou Gehrig, and Mantle batted directly before Yogi Berra. While both Ruth and Mantle drew plenty of walks, presumably pitchers rarely chose to walk them intentionally just to face Gehrig or Berra.

      I was also curious about their defensive metrics. For a fair comparison, I filtered the data to isolate only their outfield appearances. There is no value in comparing an outfield throw to an assist Ruth made on a comeback grounder while pitching. Similarly, I excluded Mantle’s infield appearances; he primarily played first base in his final two seasons and filled in briefly at other infield spots early in his career.

      Ruth played 222 more outfield games than Mantle (2,241 versus 2,019). Their total outfield putouts were nearly identical, likely because Mantle played centerfield and covered more ground. However, Ruth recorded nearly twice as many outfield assists (204 versus 117), which aligns with the arm strength expected of a former pitcher. While Ruth’s 204 assists don’t match stars such as Roberto Clemente’s 266, it remains a highly respectable number.

      Ultimately, the data shows that even if you completely ignore his pitching career, Babe Ruth built an outstanding, standalone career as both a hitter and a fielder.

              
          MANTLE  RUTH
games       2401  2503
at_bats     8102  8398
runs        1677  2174
hits        2415  2873
doubles      344   506
triples       72   136
home_runs    536   714
BA         0.298 0.342
rbi         1509  2217
sb           153   123
bb          1733  2062
so          1710  1330
games_of    2019  2241
putouts     4438  4444
assists      117   204

      Here is my R code:

library(Lahman)
library(tidyverse)

data(People)
which(People$nameLast == "Mantle")   # 11773,  playerID = mantlmi01
which(People$nameLast == "Ruth")   # 16674 , playerID = ruthba01 
data(Batting)    
data(Fielding)   # for fielding, want games in outfield POS == 'OF'

df_fielding <- Fielding %>% 
  filter(playerID %in% c('mantlmi01', 'ruthba01'), POS == 'OF') %>% 
  group_by(playerID) %>% 
  summarize(games_of = sum(G, na.rm = TRUE), putouts = sum(PO, na.rm = TRUE),
    assists = sum(A, na.rm = TRUE)
  )

df <- Batting %>% 
  filter(playerID %in% c('mantlmi01', 'ruthba01')) %>% 
   group_by(playerID) %>% 
   summarize(games = sum(G, na.rm = TRUE), at_bats = sum(AB, na.rm = TRUE),
                runs = sum(R, na.rm = TRUE), hits = sum(H, na.rm = TRUE),
                doubles = sum(X2B, na.rm = TRUE), 
                triples = sum(X3B, na.rm = TRUE), home_runs = sum(HR, na.rm = TRUE),
                BA = round(hits / at_bats,3), 
                rbi = sum(RBI, na.rm = TRUE),  sb = sum(SB, na.rm = TRUE),
                bb = sum(BB, na.rm = TRUE),  so = sum(SO, na.rm = TRUE),
                )  %>%
  merge(df_fielding, by = "playerID", all.x = TRUE) 

df_transposed <- as.data.frame(t(df))
colnames(df_transposed) <- c("MANTLE", "RUTH")
df_transposed <- df_transposed[-1, ]
df_transposed

End

Remember the “AI rat penis” incident at Frontiers in Cell & Developmental Biology?

Well, that was two and a half years ago. That’s enough water under the bridge to now check in and see how this incident affected the journal. If we look at how many papers published, we should be able to see whether the incident affected submissions.

I will say upfront that this analysis did not reveal what I thought it would. It’s well known that Frontiers journals had huge growth some years ago and I also knew that submissions had decreased in the last few years. So, I expected that this drop off was due to authors not wanting to publish in the journal after the AI figure debacle or perhaps a combination of that and the study by Mark Hanson and colleagues showing suspiciously fast turnaround times at journals, including titles from this publisher.

The plots can be quickly generated using {PubMedLagR} which is available here.

Jump to the code, or just look at the plots.

If we pull all the articles in PubMed for this journal and see how many were published per month, we get a plot like this. The dotted line shows the date of the rat penis incident.

So there was a huge boom in publications in the COVID-19 pandemic era and then a decrease. This decrease predates the rat penis incident. If anything, the number of papers published began to stabilise and perhaps even recover following the incident!

It’s a bit hard to see if the composition of articles has changed over time from the plot above. So let’s replot these data and look at the fraction of all articles that are papers versus other types. I did some spot-checking and there were several papers that are clearly Commentary, Editorial and Review types that are not properly tagged on PubMed. So take this with a pinch of salt. The articles are split 50:50 between papers and other types and this has not really changed. It could be argued that commissioned content might be unaffected whereas as directly submitted papers might drop off if authors were concerned about a journal. But there’s no evidence of that here.

This is looking at published papers. Obviously this metric lags behind the author behaviour we’re interested in. Did the rat penis incident affect submissions to the journal? We can look at the received dates of published papers:

Pretty much the same story: month-on-moth submission are steadily rising after the incident.

Obviously we only have submission data for papers that were eventually published in the journal. We don’t know what the rejection rate was over this period. It would be possible for a beleaguered journal to reduce the rejection rate to stabilise the numbers of papers it publishes. But if we assume the rejection rate is constant, then submissions follow similar dynamics.

Conclusion: the AI debacle didn’t negatively affect submissions at this journal.

Maybe potential authors were unaware of this issue (although there was plenty of publicity), or maybe they view the incident as a one-off, or perhaps they’re not concerned about science integrity issues.

It’s no secret that the biggest drivers of journal choice for authors are 1) the journal impact factor and 2) the turnaround time/hassle to publish the work. The journals with the highest impact factor and the fastest turnaround win big. A scandal like this doesn’t seem to have affected journal choice. Of course, this is just a simple analysis and there are many factors at play here.

Isn’t it ironic though that authors want to publish in a journals with high impact factors? They want their paper to be viewed as potentially highly citable because of where it is published; but if that means that their paper appeared in a journal that has a scientific integrity issue, the cry is “judge the paper on its own merits, not where it’s published”…

The code

library(PubMedLagR)
library(ggplot2)

retrieve_journal_year_records("Front Cell Dev Biol", 2015:2026, batch_size = 200, papers_only = FALSE)
all <- pubmed_xmls_to_df(clean = FALSE)
# remove duplicate rows
all <- all[!duplicated(all),]
# clear Data/ and then
retrieve_journal_year_records("Front Cell Dev Biol", 2015:2026, batch_size = 200)
pprs <- pubmed_xmls_to_df()
# classify the papers in all
all$paper <- ifelse(all$pmid %in% unique(pprs$pmid),"paper", "other")
# convert publication date to the first of the month
all$year_month <- as.Date(paste0(substr(all$pubdate,1,7),"-01"))

# frontiers colours for fun
x <- c("242 130 91", "24 157 196")
fcolours <- sapply(strsplit(x, " "), function(x)
  rgb(x[1], x[2], x[3], maxColorValue=255))

# ggplot of number of papers per month-year
p1 <- ggplot(all, aes(x = year_month, fill = paper)) +
  geom_bar() +
  scale_fill_manual(values = fcolours) +
  geom_vline(xintercept = as.POSIXct(as.Date("2024-02-14")), linetype=2) +
  labs(x = "", y = "Papers per month") +
  theme_classic()
p1

# stacked plot version
p2 <- ggplot(all, aes(x = year_month, fill = paper)) +
  geom_bar(position = "fill") +
  scale_fill_manual(values = fcolours) +
  geom_vline(xintercept = as.POSIXct(as.Date("2024-02-14")), linetype=2) +
  labs(x = "", y = "Composition") +
  theme_classic()
p2

# what about submissions?
all$sub_year_month <- as.Date(paste0(substr(all$recdate,1,7),"-01"))
p3 <- ggplot(all, aes(x = sub_year_month, fill = paper)) +
  geom_bar() +
  scale_fill_manual(values = fcolours) +
  geom_vline(xintercept = as.POSIXct(as.Date("2024-02-14")), linetype=2) +
  labs(x = "", y = "Papers per month") +
  theme_classic()
p3

I added some styling to the plots to post them here, but the code will give the same three plots as in the post.

—

The post title comes from “New Frontier” by Donald Fagen from his 1982 “The Nightfly” LP.

controller is an R package for working with controlled vocabularies. It’s first release (v0.1.0) is now available now on CRAN.

The package addresses something I find myself doing very often in analysis code: tidying messy and inconsistent terminologies. For smaller datasets, dplyr::recode() is okay for this, but writing the mapping out as an R function call gets tedious fast when dealing with a long list of terms. It becomes very tedious when you have variants distinguished only by things like capitalisation (OxA- vs. oxa-), word boundaries (Çatalhöyük vs. Çatal Höyük) or character encoding (ʿAin Ghazal vs. ʽAyn Ghazal).

controller instead defines preferred terms and their variants in a data frame. Its control() verb is the equivalent of dplyr::recode() but using this thesaurus and with a few extra bells and whistles for fuzzy matching and reporting what was (and wasn’t) changed:

library(controller)
data("colour_thesaurus")

shades <- c("daffodil", "purple", "magenta", "azure", "navy", "violet")
control(shades, colour_thesaurus)
#> Replaced values:
#> ℹ daffodil → yellow
#> ℹ azure → blue
#> ℹ navy → blue
#> ℹ violet → purple
#> Warning: Some values of `x` were not matched in `thesaurus`:
#> ✖ magenta

Fuzzy matching means we don’t need to exhaustively list those variants from things like differences in case, word boundaries, or character encoding:

control_ci(toupper(shades), colour_thesaurus)
#> Replaced values:
#> ℹ DAFFODIL → yellow
#> ℹ PURPLE → purple
#> ℹ AZURE → blue
#> ℹ NAVY → blue
#> ℹ VIOLET → purple
#> [1] "yellow"  "purple"  "MAGENTA" "blue"    "blue"    "purple"
#> Warning message:
#> Some values of `x` were not matched in `thesaurus`:
#> ✖ MAGENTA

This package has been hanging around for a while! It started off as a helper function I used for cleaning up site names from prehistoric sites in Southwest Asia. The basic idea was inspired by similar functions that used to exist in c14bazAAR for cleaning sample metadata for radiocarbon date, that I thought were quite neat. So when the maintainers of that package decided to deprecate those, I took over the thesauri as part of c14 and spun the supporting functions off into controller as a standalone package. Then over the years it acquired some more functionality for working with controlled vocabularies (a surprising gap in the R ecosystem), like reading heritage vocabularies in Historic England’s FISH format. Five years later, I am finally getting around to releasing it on CRAN because I need to release c14 on CRAN, because that’s used in analyses I’m now publishing. It’s the research software engineering of changing a lightbulb, basically.

The first release of controller includes:

• control(), control_ci(), and control_fuzzy() for recoding values
• control_names(), control_names_ci(), and control_names_fuzzy() for recoding names
• control_matches() for inspecting how matches were made
• read_fish() for reading vocabularies in Historic England’s FISH format
• colour_thesaurus, an example dataset

You can install it from CRAN:

install.packages("controller")

Or the development version from GitHub:

remotes::install_github("joeroe/controller")

You can find the full documentation at https://controller.joeroe.io.

muttest so far only printed to the console which mutants were killed and which survived.

That’s the minimal setup you need to act on mutation testing findings, but it’s far from optimal. Clear the terminal, close the IDE, and the results are gone unless you saved the output.

muttest@0.3.0 fixes that with an HTML report.

1. Run your tests with new JSON reporter.
2. Render the report.
3. Open it in your browser, and every mutant is laid over your source code.

Click the ones that survived and see exactly what slipped past your suite. That list of survivors is your to-do list: every one is a candidate test. You decide which are worth writing and which are just noise.

Read the full changelog here.

New to mutation testing?

Coverage tells you which lines ran. It says nothing about whether your tests would notice if those lines were wrong. You can delete every assertion in your suite, run covr, and still see 100%.

Mutation testing asks the harder question: if this code were subtly broken, would a test fail? It makes small changes to your source (> becomes >=, TRUE becomes FALSE) and reruns your tests against each version. Each changed version is a mutant. If a test fails, the mutant is killed. If every test still passes, the mutant survived, and you’ve found a gap.

A surviving mutant isn’t a bug to fix. It’s a missing test that could hide an expensive bug later.

If you want to dive deeper, the 0.2.0 post walks through it with an example.

The report

Run your mutation tests with JSONMutationReporter, call report(), and you get a self-contained HTML file.

See the full example report here.

Here’s what the report does that the console can’t do:

• Your source code annotated with mutants. Every file is shown with its actual code, and each mutated line is marked inline. You read the survivor in context instead of reconstructing where a given line was.
• Per-file scores. The overall number is a starting point, not a diagnosis. The report breaks the score down by file, so you can see your parser sitting at 95% while the validation module is at 40%. That tells you which file to open first.
• Per-line mutant diffs. Click a mutant and see exactly what changed: the original line, the mutated line, and whether your tests killed it. When something survived, the diff tells you what kind of a potential bug slipped through.
• Filter by status. The report opens focused on survivors, since the killed mutants already did their job. You can bring the others back when you want the full picture, but the default is the view you actually came for.
• Keyboard navigation. Step through mutants without scrolling through the whole file.

That’s the workflow – filter to survivors, walk the list top to bottom, and each one is a decision: write the test, or call it noise and move on.

A checklist you can work through.

JSON output

The new JSONMutationReporter writes results to a file that conforms to the mutation-testing-elements schema. That schema isn’t something I invented. It’s the format shared by StrykerJS, Stryker.NET, Stryker for Scala, and others. If other teams in your organization already use Stryker and its reporting tools, you can plug muttest results straight into the same dashboards.

library(muttest)

plan <- muttest_plan(
  source_files = "R/shipping.R",
  mutators = comparison_operators()
)

muttest(plan, reporter = JSONMutationReporter$new())

report()

JSONMutationReporter writes to muttest.json by default, and report() reads that same file and writes muttest.html next to it, so with the defaults there’s nothing to wire up.

If you don’t like the built-in report, the JSON file is there to build your own. I plan to keep this schema stable so don’t worry about your report breaking with every next release.

Watch it run, save it for later

You usually want both: a live progress table while the run happens, and a saved file to open afterward. Before, you picked one reporter and got one behavior.

MultiReporter runs several reporters at once. Feed it a ProgressMutationReporter and a JSONMutationReporter together. The console fills up as mutants are tested, and the JSON file is waiting when it’s done.

muttest(
  plan,
  reporter = MultiReporter$new(
    ProgressMutationReporter$new(),
    JSONMutationReporter$new()
  )
)

It’s exactly the same idea as testthat’s MultiReporter.

Scoring changes worth knowing about

Two changes in this release will move your score:

• Errored mutants now count as killed. Some mutations don’t just fail an expectation; they make the code throw before a test can even assert anything. A mutation that produces an error is still a mutation your suite detected, so 0.3.0 scores it as killed. This matches PIT, Stryker, and mutmut. Previously muttest counted these as survived, which was too harsh: the test did catch the change, it just caught it by blowing up.
• Uncovered code gets its own category. Under FileTestStrategy, a mutant in a source file with no matching test file used to look identical to a mutant your tests ran and missed. Those are different problems. A survived mutant means your tests are weak; an uncovered mutant means there are no tests at all. 0.3.0 reports uncovered mutants as a separate no coverage category and leaves them out of the score, so untested code stops masquerading as escaped mutants and dragging your number down for the wrong reason.

If your score shifts after upgrading, this is why. It’s more honest this way.

Tell me about your experience of using muttest

I hope this release makes it easier to analyze results.

muttest is still young and the interface may shift. If the report is missing something you’d want to see, if a mutation you care about has no way to be expressed, or if the new scoring surprises you in a way that feels wrong, open an issue on GitHub.

Feature requests are just as welcome as bug reports.

We all were presented with this problem – water jug problem, which – in the time of Football world cup 2026 – can be translated to any liquid. *hint hint* But the riddle is as logical as mathematical. With mathematics finding the greatest common divisor. In general, it can be used with state search or Depth First Search (DFS).

With DFS we can solve this with Bezout’s identity; which in general is a theorem which relates two arbitraty integers with their greatests common divisor; and used in algebraic language, finding common zeros of n-polznomials in n-indeterminates. So the common zeros equals the product of the degrees of the polynomials.

And now that the imagine splitting the 16L beer Giraffe in two 8L Giraffes but you are only using 11L and 7L empty giraffes

This is the proof that with beer, Algebra is more fun And because it is fun, we can also find greated common divisors using Breadth-First Search (BFS). And here is the code:

solve_jugs <- function(caps = c(16, 11, 7), start = c(16, 0, 0), goal = c(8, 8, 0)) {
  # BFS over all (a, b, c) states
  # Each state is a named integer vector of water = amount is each jug or ?????
  
  queue   <- list(list(state = start, path = list(start)))
  visited <- list()
  key     <- function(s) paste(s, collapse = "-")
  
  while (length(queue) > 0) {
    node  <- queue[[1]]
    queue <- queue[-1]
    s     <- node$state
    
    if (isTRUE(all(s == goal))) return(node$path)
    if (!is.null(visited[[key(s)]])) next
    visited[[key(s)]] <- TRUE
    
    n <- length(s)
    for (from in 1:n) {
      for (to in 1:n) {
        if (from == to || s[from] == 0 || s[to] == caps[to]) next
        
        pour    <- min(s[from], caps[to] - s[to])
        new_s   <- s
        new_s[from] <- s[from] - pour
        new_s[to]   <- s[to]   + pour
        
        if (is.null(visited[[key(new_s)]])) {
          queue <- c(queue, list(list(
            state = new_s,
            path  = c(node$path, list(new_s))
          )))
        }
      }
    }
  }
  NULL  # no solution; add message or smht :)
}

solution <- solve_jugs()

 
for (step in solution) {
  cat(sprintf("  %-4d  %-4d  %-4d\n", step[1], step[2], step[3]))
}
  

And the final solution will reveal the steps and actions:

  16    0     0   
  5     11    0   
  5     4     7   
  12    4     0   
  12    0     4   
  1     11    4   
  1     8     7   
  8     8     0   
  

Similar steps are presented on the animation above.

As always, the complete code is available on GitHub in  Useless_R_function repository. And code to the animation is here (same Github repository).

Check the repository for future updates!

Stay healthy, hydrated and happy R-coding!

At useR! 2026 in Warsaw, I contributed to a workshop on Geocomputation with R (July 6, 2026) and gave a keynote entitled A world still to be mapped: reflections on geocomputation in R (July 8, 2026). Both were opportunities to consider what the R spatial ecosystem already makes possible, its strengths, limitations, and what remains to be done.

Keynote slides: https://jakubnowosad.com/user2026/

Workshop materials: https://github.com/geocompx/user26

@online{nowosad2026,
  author = {Nowosad, Jakub},
  title = {A World Still to Be Mapped: Reflections on Geocomputation in
    {R:} Takeaways from the Talk and Workshop at {UseR!} 2026},
  date = {2026-07-21},
  url = {https://jakubnowosad.com/posts/2026-07-21-user/},
  langid = {en}
}

TLDR: The first post established that the WGI and WJP rule of law measures agree on how they rank countries. This follow-up tests the stricter question of whether the WGI can predict the WJP’s levels, rather than just the relative ranking / relative position of countries. It can: forecasting errors appear to be minimal and don’t vary significantly over time.

This is second post in a series examining relationships between the rule of law and economic and social outcomes.

In the previous post in the series I tested whether the World Governance Indicator’s (WGI) rule of law index (RoL) was a worthy substitute for the index produced by the World Justice Project when needing to take advantage of the WGI’s wider and longer coverage.

When I asked Claude for feedback it rudely provided it, arguing that my analysis was fine for confirming general agreement as to how the two measures ranked countries, but additional analysis was needed to test whether the WGI could reliably be used to predict RoL levels.

My first response was to tell Claude that nobody should be using composite indices in this way in the first place as it implied a level of precision they don’t have, but I decided it would be a useful addendum to the first post that I could post during my holiday in the Northern Territory.

Project Setup and Data

Once again, data used in this post can be downloaded here for the WGI and here for the WJP’s RoL index. These were current as of July 2025.

#load the packages we'll probably need
library(tidyverse)
library(readxl)
library(janitor)
library(countrycode)

#import WGI data
dta_wgi_2025<-read_excel("./Data/wgidataset_with_sourcedata-2025.xlsx",
                         sheet="rl") |> 
  clean_names()

#import World Justice Project RoL data
dta_wjp_rol<-read_excel('./Data/2025_wjp_rule_of_law_index_HISTORICAL_DATA_FILE.xlsx', sheet="Historical Data")|> 
  clean_names()

Data Cleaning

Data cleaning is identical to the last post and revolves around making sure that countries names and regions are consistently applied across the two indices.

#standardize column names
dta_wgi_2025<-dta_wgi_2025 |> 
  rename(country=economy_name,
         iso3c=economy_code,
         wgi_rol=governance_estimate_approx_2_5_to_2_5)

dta_wjp_rol<-dta_wjp_rol |> 
  rename(iso3c=country_code) |> 
  rename_with(~ str_replace(., "^x", "factor_"), starts_with("x"))

#change wjp's year variable to YYYY format and convert to numeric 
#(adopts the first 4 digit year when in YYYY-YYYY format)
dta_wjp_rol <- dta_wjp_rol |> 
  mutate(year = as.numeric(str_sub(year, 1, 4)))

#cold-heartedly drop columns I'm not interested in 
dta_wgi_2025<-dta_wgi_2025 |> 
  select(iso3c, income_classification, year,wgi_rol)
#drop country and region name labels so these can be standardized 
dta_wjp_rol<-dta_wjp_rol |> 
  select(-country_year,-country,-region) |> 
  rename(wjp_rol=wjp_rule_of_law_index_overall_score)


#merge dataframes
dta_rol_unified<-left_join(dta_wgi_2025,
                           dta_wjp_rol,
                           by = join_by(year, iso3c), 
                           keep=FALSE) 

#add standardized and region names country names

#define country code assignments for legacy / ambigious codes
#(Note: matches devised by Claude) 
ref_iso3c_custom_names <- c(ADO = "Andorra",
                            ANT = "Netherlands Antilles",
                            PRI = "Puerto Rico",
                            REU = "Réunion",
                            XKX = "Kosovo")

ref_iso3c_custom_regions <- c(ADO = "Europe & Central Asia",
                              ANT = "Latin America & Caribbean",
                              PRI = "Latin America & Caribbean",
                              REU = "Sub-Saharan Africa",
                              XKX = "Europe & Central Asia")

#assign country names and regions:
dta_rol_unified<-dta_rol_unified |> 
  mutate(country_name=countrycode(iso3c, 
                                  origin='iso3c',
                                  destination = 'country.name.en',    
                                  custom_match = ref_iso3c_custom_names),
         region=countrycode(iso3c, 
                            origin='iso3c',
                            destination = 'region',    
                            custom_match =ref_iso3c_custom_regions))

Testing Cardinal Substitution

The code below basically asks whether it’s possible to estimate the level of WJP’s RoL index using the WGI’s RoL measure using a simple linear model. At the outset, the model’s fit is quite high at 97 percent, but we could have guessed this from the cross-country scatter we produced in our last post.

# Cardinal substitution: does a WGI->WJP mapping reproduce WJP levels? -----
sum_rol_country <- dta_rol_unified |>
  filter(!is.na(wgi_rol), !is.na(wjp_rol)) |>
  summarise(wgi_rol_mean = mean(wgi_rol),
            wjp_rol_mean = mean(wjp_rol),
            .by = country_name)

mod_rol_calib <- lm(wjp_rol_mean ~ wgi_rol_mean, data = sum_rol_country)

#output the R squared
summary(mod_rol_calib)$r.squared |> round(2)

The second chunk calculates the Leave-One-Out Cross-Validation Root Mean Squared Error (LOO RMSE), which estimates prediction error for each country when it’s held out of the sample: a smaller LOO RMSE implies smaller forecast errors when predicting outside the sample. This is then scaled by the standard deviation of the country-level WJP means to make it easier to interpretat.

The results are encouraging: the LOO RMSE indicates that using the WGI to estimate WJP levels carries a typical error of less than three percent, which is roughly 19% of the cross-country spread in WJP scores. Suggesting the WGI can provide a reasonable proxy when the levels are important too.

# Leave-One-Out Cross-Validation Root Mean Squared Error (LOO RMSE)
rlt_rol_loo_rmse   <- sqrt(mean((residuals(mod_rol_calib) /
                                   (1 - hatvalues(mod_rol_calib)))^2))

#divide by standard deviation of WJP index means to put the figure in the context fo the data
rlt_rol_rmse_ratio <- rlt_rol_loo_rmse / sd(sum_rol_country$wjp_rol_mean)

Systematic Bias

Of course, because the LOO RMSE measures is calculated across the entire dataset, it might hide systematic bias, such as if the forecast accuracy varies depending on the WGI score.

This is tested in the plot below, by comparing the prediction error (residual) with the predicted value for the WJP’s index (using the WGI RoL measure). If there no systematic bias, the residuals should remain relatively stable across the predicted WJP values.

Alas, this doesn’t appear to be the case, with the model under-predicting the WJP’s RoL index at both ends. Although this bias is relatively mild at 0.05 (in WJP index units), it does suggest the WGI is likely to be a poorer substitute for the WJP’s index for countries at either end of the RoL spectrum.

dta_plt_rol_calib <- sum_rol_country |>
  mutate(wjp_pred = fitted(mod_rol_calib),
         resid    = wjp_rol_mean - wjp_pred)

plt_rol_calib <- ggplot(dta_plt_rol_calib, aes(wjp_pred, resid)) +
  geom_hline(yintercept = 0, linewidth = 0.3, col = "grey70") +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, col = "grey40", linewidth = 0.5) +
  labs(x = "Predicted WJP (from WGI)", y = "Residual (WJP units)") +
  theme_classic()

plt_rol_calib

However, being an applied economist means I know a simple trick: add non-linear terms to everything and see if the smile disappears.

And it does, the smile is now a grimace: as is shown in the plot below adding a non-linear term seems to take care of the bias at either end. Notice the average errors no longer increase at either end, making the prediction relatively constant (and small) across the RoL spectrum. Suggesting the WGI can be trusted to predict the WJP index at either end of the spectrum provided we account for the non-linearity.

#reestimate the model with a non-linear WGI term 
mod_rol_calib_non_linear <- lm(wjp_rol_mean ~ wgi_rol_mean + I(wgi_rol_mean^2),
                               data = sum_rol_country)


dta_plt_rol_calib_nl <- sum_rol_country |>
  mutate(wjp_pred = fitted(mod_rol_calib_non_linear),
         resid    = wjp_rol_mean - wjp_pred)

plt_rol_calib_nl <- ggplot(dta_plt_rol_calib_nl, aes(wjp_pred, resid)) +
  geom_hline(yintercept = 0, linewidth = 0.3, col = "grey70") +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "loess", se = FALSE, col = "grey40", linewidth = 0.5) +
  labs(x = "Predicted WJP (from WGI, quadratic)", y = "Residual (WJP units)") +
  theme_classic()

plt_rol_calib_nl

But does it drift?

The final question Claude suggested I ask was whether the WGI’s ability to predict the WJP’s index changes over time. The code below tests this by comparing a simple linear prediction model with one where the influence of the WGI changes over time.

The anova() function then compares results of both models to test whether allowing the WGI’s influence to change over time improves the model’s forecasting accuracy. If it does (indicated by a significant p-value), the WGI’s ability to predict the WJP’s index might not be stable over time.

The result: a p value of 9%, which points to the time invariant model being good enough for predicting the WJP’s index.

#estimate a simple linear model with WGI's influence kept constant over time
mod_rol_calib_flat <- lm(wjp_rol ~ wgi_rol, data = dta_rol_unified)

#estimate a model where WGI's influence varies over time
mod_rol_calib_yr   <- lm(wjp_rol ~ wgi_rol * factor(year), data = dta_rol_unified)

#compare the prediction capacity of both models 
anova(mod_rol_calib_flat, mod_rol_calib_yr)

Look what Claude made me do…

So there you have it: the WGI can serve as a useful substitute for the WJP’s rule of law index whether you’re interested in RoL rankings or levels, provided you account for the non-linearity at either end of the RoL scale.

How AI was used for this post: As this additional analysis was Claude’s suggestion and I was hours away from a holiday I leaned heavily on Claude for the first draft of the code. Much of the text is mine, but I did ask AI for suggestions to make sure I was explaining the analysis accurately.

The post Level With Me: Can the WGI Predict the WJP’s Rule of Law Scores? appeared first on Giles.

Usually, the first step in every nonlinear regression analysis is to select the function that best describes the phenomenon under study. The next step is to fit this function to the observed data, possibly by using some sort of nonlinear least squares algorithm.

We have already devoted a post to the first task, which you can find at this link.

As for the second step, the main problem is that nonlinear least squares algorithms are iterative, in the sense that they start from some initial guesses for the model parameters, which are continuously improved until the least squares solution is approximately reached. Quite often, providing such initial guesses for all model parameters becomes a problem: if our guesses are not close enough to the least squares estimates, the algorithm may stall and fail to converge. Or, even worse, it may converge to the wrong solution. How do we obtain good initial guesses for the model parameters? This is not easily accomplished, especially for students and practitioners. This is where self-starters come in handy.

Self-starter functions can automatically calculate initial values for any given dataset and, therefore, they can make nonlinear regression almost as straightforward as linear regression. From a teaching perspective, this means that the transition from linear to nonlinear models is immediate and hassle-free.

In another post, at this link, I explained how self-starters can be built for both the nls() function in the ‘stats’ package and the drm() function in the ‘drc’ package (Ritz et al., 2019). In this post, I would like to provide an overview of the self-starting functions that already exist in R, either in the ‘stats’, ‘drc’, or ‘statforbiology’ packages. I do not aim for completeness here, as other packages also contain self-starters, such as the ‘nlraa’ package (Miguez, 2025), with which I am not sufficiently familiar. The exemplary datasets are included in the ‘statforbiology’ package and they all come from real field or greenhouse experiments; further information and citations to the original works can be found in my book (Onofri, 2026), or in the additional meterial at this link.

Functions and curve shapes

As in the other post, I have classified linear/nonlinear regression functions according to the shape they show when they are plotted in an x-y graph, following the approach taken in Ratkowsky (1990):

• Polynomials
  1. Straight line function
  2. Quadratic polynomial function
• Concave/Convex curves (no inflection)
  1. Exponential function
  2. Asymptotic function / Negative exponential function
  3. Power function
  4. Logarithmic function
  5. Rectangular hyperbola
• Sigmoidal curves
  1. Logistic function
  2. Gompertz function
  3. Modified Gompertz function
  4. Log-logistic function
  5. Weibull (type-1) function
  6. Weibull (type-2) function
• Curves with maxima/minima
  1. Peaked sigmoidal function
  2. Bragg function
  3. Lorentz function
  4. Beta function

To make navigation easier, you can first inspect Figure 1 to identify the type of curve of interest and then use the links above to jump directly to the corresponding section.

First of all, we need to install (if necessary) and load these packages, by using the code below.

# installing package, if not yet available
# install.packages("statforbiology")

# loading package
library(statforbiology)

# Example of fitting four straight lines at once
library(statforbiology)
dataset <- getAgroData("metamitron")
model <- drm(log(Conc) ~ Time, fct = DRC.linear(),
             curveid = Herbicide, data = dataset)
summary(model)

Model fitted: Straight line (2 parms)

Parameter estimates:

        Estimate Std. Error t-value   p-value    
a:M    4.5159111  0.0577603  78.184 < 2.2e-16 ***
a:MP   4.6304853  0.0577603  80.167 < 2.2e-16 ***
a:MC   4.5401707  0.0577603  78.604 < 2.2e-16 ***
a:MPC  4.7546974  0.0577603  82.318 < 2.2e-16 ***
b:M   -0.0386165  0.0015585 -24.777 < 2.2e-16 ***
b:MP  -0.0318908  0.0015585 -20.462 < 2.2e-16 ***
b:MC  -0.0231879  0.0015585 -14.878 < 2.2e-16 ***
b:MPC -0.0234014  0.0015585 -15.015 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error:

 0.1687429 (88 degrees of freedom)

library(statforbiology)
dataset <- getAgroData("waterProductivity")

# nls fit
model <- nls(WP ~ NLS.poly2(Regime, a, b, c),
             data = dataset)

#drc fit
model <- drm(WP ~ Regime, fct = DRC.poly2(),
             data = dataset)

Model fitted: Second Order Polynomial (3 parms)

Parameter estimates:

                 Estimate  Std. Error  t-value   p-value    
a:(Intercept)  1.9187e+00  2.3171e-01   8.2806 1.679e-05 ***
b:(Intercept)  4.2977e-01  8.4554e-03  50.8287 2.217e-12 ***
c:(Intercept) -3.3618e-03  6.6586e-05 -50.4884 2.355e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error:

 0.1441625 (9 degrees of freedom)

library(statforbiology)
data(degradation)

# nls fit
model <- nls(Conc ~ NLS.expoDecay(Time, a, k),
             data = degradation)

# drm fit
model <- drm(Conc ~ Time, fct = DRC.expoDecay(),
             data = degradation)
summary(model)

Model fitted: Exponential Decay Model (2 parms)

Parameter estimates:

                 Estimate Std. Error t-value   p-value    
C0:(Intercept) 99.6349312  1.4646680  68.026 < 2.2e-16 ***
k:(Intercept)   0.0670391  0.0019089  35.120 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error:

 2.621386 (22 degrees of freedom)

Time <- c(1, 3, 5, 7, 9, 11, 13, 20)
Weight <- c(8.22, 14.0, 17.2, 16.9, 19.2, 19.6, 19.4, 19.6)

# nls fit
model <- nls(Weight ~ NLS.asymReg(Time, b, m, a) )

# drm fit
model <- drm(Weight ~ Time, fct = DRC.asymReg(names = c("b", "m", "a")))
summary(model)

Model fitted: Asymptotic Regression Model (3 parms)

Parameter estimates:

               Estimate Std. Error t-value   p-value    
b:(Intercept)  3.756123   1.367835  2.7460  0.040500 *  
m:(Intercept)  0.337078   0.052704  6.3957  0.001385 ** 
a:(Intercept) 19.629817   0.420557 46.6758 8.525e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error:

 0.639656 (5 degrees of freedom)

library(statforbiology)
speciesArea <- getAgroData("speciesArea")

#nls fit
model <- nls(numSpecies ~ NLS.powerCurve(Area, a, b),
             data = speciesArea)

# drm fit
model <- drm(numSpecies ~ Area, fct = DRC.powerCurve(),
             data = speciesArea)
summary(model)

Model fitted: Power curve (Freundlich equation) (2 parms)

Parameter estimates:

              Estimate Std. Error t-value   p-value    
a:(Intercept) 4.348404   0.337197  12.896 3.917e-06 ***
b:(Intercept) 0.329770   0.016723  19.719 2.155e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error:

 0.9588598 (7 degrees of freedom)

library(statforbiology)
speciesArea <- getAgroData("speciesArea")

#nls fit
model2 <- nls(numSpecies ~ NLS.logCurve(Area, a, b),
             data = speciesArea)

# drm fit
model2 <- drm(numSpecies ~ Area, fct = DRC.logCurve(),
             data = speciesArea)
summary(model2)

Model fitted: Linear regression on log-transformed x (2 parms)

Parameter estimates:

              Estimate Std. Error t-value   p-value    
a:(Intercept)  1.51111    1.17401  1.2871     0.239    
b:(Intercept)  4.06359    0.35575 11.4224 8.847e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error:

 1.910082 (7 degrees of freedom)

       df      AIC
model   3 28.52288
model2  3 40.92769

library(statforbiology)
dataset <- getAgroData("Ammi94_YL")

#drm fit
model <- drm(YL ~ Density, fct = DRC.YL(), data = dataset)
summary(model)

Model fitted: Yield-Loss function (Cousens, 1985) (2 parms)

Parameter estimates:

               Estimate Std. Error t-value p-value  
i:(Intercept) 0.0118359  0.0045637  2.5935 0.02228 *
A:(Intercept) 0.4841796  0.1863969  2.5976 0.02211 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error:

 0.05453432 (13 degrees of freedom)

Sigmoidal curves

In the previous post (at this link) we examined six different sigmoidal shapes based on , to be selected according to the expected symmetry of the response. The equations are the logistic, Gompertz, and modified Gompertz functions, with the following three equations, respectively:

The same curves can be based on the logarithm of , to obtain the corresponding log-logistic, Type-1 Weibull, and Type-2 Weibull equations, which are, respectively:

The meaning of the model parameters and details about the shapes can be found in another post at this link.

For all these four-parameter curves, it is possible to set to obtain the three-parameter curves (17a), (18a), (19a), (20a), (21a), and (22a). Likewise, we can constrain both and to obtain the two-parameter curves (17b), (18b), (19b), (20b), (21b), and (22b). The collection of self-starters is huge; therefore, I will provide a summary table.

List of sigmoidal functions, together with self-starters for fitting with drm() and nls(). (1) These self-starters are included in the ‘drc’ package and the corresponding curves have an increasing shape when and a decreasing shape when . For all other self-starters, which are available in the ‘statforbiology’ package, the reverse is true (which is the most common definition, apart from the case of herbicide bioassays){#tbl-models}

Name	Eq. No.	for drm()	for nls()
Logistic (4-parameters)	(17)	L.4() (1)	NLS.L4()
Logistic (3-parameters)	(17a)	L.3() (1)	NLS.L3()
Logistic (2-parameters)	(17b)	L.2() (1)	NLS.L2()
Gompertz (4-parameters)	(18)	G.4() (1)	NLS.G4()
Gompertz (3-parameters)	(18a)	G.3() (1)	NLS.G3()
Gompertz (2-parameters)	(18b)	G.2() (1)	NLS.G2()
mod-Gompertz (4-parameters)	(19)	DRC.E4()	NLS.E4()
mod-Gompertz (3-parameters)	(19a)	DRC.E3()	NLS.E3()
mod-Gompertz (2-parameters)	(19b)	DRC.E2()	NLS.E2()
Log-logistic (4-parameters)	(20)	LL.4() (1)	NLS.LL4()
Log-logistic (3-parameters)	(20a)	LL.3() (1)	NLS.LL3()
Log-logistic (2-parameters)	(20b)	LL.2() (1)	NLS.LL2()
Type-1 Weibull (4-parameters)	(21)	W1.4() (1)	NLS.W1.4()
Type-1 Weibull (3-parameters)	(21a)	W1.3() (1)	NLS.W1.3()
Type-1 Weibull (2-parameters)	(21b)	W1.2() (1)	NLS.W1.2()
Type-2 Weibull (4-parameters)	(22)	W2.4() (1)	NLS.W2.4()
Type-2 Weibull (3-parameters)	(22a)	W2.3() (1)	NLS.W2.3()
Type-2 Weibull (2-parameters)	(22b)	W2.2() (1)	NLS.W2.2()

Example 8

The dataset ‘beetGrowth’ in the ‘statforbiology’ package reports the results of an experiment in which the growth of sugar beet was evaluated under either weed-free or weed-infested conditions. In the following example, we consider the time course of crop weight under infested conditions and fit the logistic, Gompertz, and modified Gompertz functions to the observed data. We see that, based on the AIC, the three fits are equivalent in practice. We fit the same models with ‘drm()’ to show that the two fitting functions give the same results, apart from the sign of the parameter, which is reversed in the L.3() function.

library(statforbiology)
beetGrowth <- getAgroData("beetGrowth")

# nls fit
model.1 <- nls(Infested ~ NLS.L3(DAE, b, d, e), data = beetGrowth)
model.2 <- nls(Infested ~ NLS.G3(DAE, b, d, e), data = beetGrowth)
model.3 <- nls(Infested ~ NLS.E3(DAE, b, d, e), data = beetGrowth)
AIC(model.1, model.2, model.3)

        df      AIC
model.1  4 250.2906
model.2  4 250.7596
model.3  4 250.8456

summary(model.1)

Formula: Infested ~ NLS.L3(DAE, b, d, e)

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
b 1.188e-01  1.874e-02   6.338 1.33e-05 ***
d 2.512e+03  1.277e+02  19.671 4.01e-12 ***
e 5.803e+01  1.823e+00  31.836 3.46e-15 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 222 on 15 degrees of freedom

Number of iterations to convergence: 6 
Achieved convergence tolerance: 4.299e-06

# drm fit
model <- drm(Infested ~ DAE, fct = L.3(), data = beetGrowth)
summary(model)

Model fitted: Logistic (ED50 as parameter) with lower limit fixed at 0 (3 parms)

Parameter estimates:

                Estimate Std. Error t-value   p-value    
b:(Intercept)   -0.11876    0.01832  -6.483 1.033e-05 ***
d:(Intercept) 2511.95664  127.95600  19.631 4.131e-12 ***
e:(Intercept)   58.03099    1.83458  31.632 3.788e-15 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error:

 221.9644 (15 degrees of freedom)

Example 9

The dataset ‘brassica’ in the ‘statforbiology’ package reports the results of a short-term pot bioassay to determine the phytotoxicity of imazethapyr in hydroponic solution (Onofri, 1994). The response variable is the weight of the test plant (Brassica rapa), as affected by the herbicide concentration in the nutrient solution, which requires a decreasing sigmoidal curve. In the code below, we fit both a 4-parameter log-logistic and a 3-parameter log-logistic (with ) function and graphically compare the two fits.

library(statforbiology)
brassica <- getAgroData("brassica")
model.1 <- drm(FW ~ Dose, fct = LL.4(), data = brassica)
model.2 <- drm(FW ~ Dose, fct = LL.3(), data = brassica)

library(statforbiology)
dataset <- getAgroData("SOLNI_hormesis")
model <- drm(FW ~ Dose, fct = BC.5(),
            data = dataset)
summary(model)

Model fitted: Brain-Cousens (hormesis) (5 parms)

Parameter estimates:

              Estimate Std. Error t-value   p-value    
b:(Intercept)  3.57034    1.25309  2.8492 0.0082868 ** 
c:(Intercept)  3.76335    0.85964  4.3778 0.0001619 ***
d:(Intercept) 14.14991    1.36933 10.3334 6.992e-11 ***
e:(Intercept)  2.85405    0.76931  3.7099 0.0009485 ***
f:(Intercept)  5.26894    3.38238  1.5578 0.1309350    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error:

 2.97547 (27 degrees of freedom)

library(statforbiology)
dataset <- getAgroData("Hordeum_GR")
model <- drm(GR50 ~ Temp, fct = DRC.bragg.4(),
            data = dataset)
summary(model)

Model fitted: Bragg equation with four parameters

Parameter estimates:

                Estimate Std. Error t-value   p-value    
b:(Intercept)  0.0052849  0.0036063  1.4654 0.2804269    
c:(Intercept)  0.0039006  0.0131491  0.2966 0.7947090    
d:(Intercept)  0.0290870  0.0010552 27.5646 0.0013135 ** 
e:(Intercept) 18.0013282  0.5220972 34.4789 0.0008401 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error:

 0.001544712 (2 degrees of freedom)